{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "49e3f9ef",
   "metadata": {},
   "source": [
    "# Stability of Runge-Kutta Method\n",
    "**강좌**: *수치해석*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e9d3762",
   "metadata": {},
   "source": [
    "## 안정성\n",
    "시간 간격에 따라 정확도 뿐 아니라 알고리즘이 발산할 수 있다. 이는 수치 안정성과 관계된다.\n",
    "\n",
    "안정성을 분석하기 위해 Model problem을 생각한다.\n",
    "\n",
    "$$\n",
    "\\frac{dy}{dt}=f(y,t) = f(y_0, t_0) + (t-t_0)f_t + (y - y_0) f_y + ... = \\lambda y + Others\n",
    "$$\n",
    "$$\n",
    "\\frac{dy}{dt}= \\lambda y, ~~~~\\lambda = \\lambda_R + i \\lambda_i\n",
    "$$\n",
    "\n",
    "$\\lambda_R < 0$ 인 경우 이 미분방정식은 해가 제한되어 (bounded) 있다.\n",
    "\n",
    "### Sympy\n",
    "[sympy](https://docs.sympy.org/latest/index.html) 는 Symbolic Math 프로그램으로 수식 연산을 지원한다.\n",
    "\n",
    "이를 이용하면 수식 전개에 활용할 수 있다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "851c075b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.dpi'] = 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d551c620",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "\n",
    "# Symbolic Math\n",
    "lam = sy.Symbol('lambda')\n",
    "h = sy.Symbol('h')\n",
    "t = sy.Symbol('t')\n",
    "y = sy.Symbol('y')\n",
    "\n",
    "\n",
    "# Define model function y'= lambda y\n",
    "def ff(t, y):\n",
    "    return lam*y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e59eebb",
   "metadata": {},
   "source": [
    "### Explicit Euler Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d1dea78b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h*lambda + 1\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Stability region of Explicit Euler')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzcAAAKOCAYAAACbY+fTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAABKH0lEQVR4nO3deZyVZd0/8O+wr7JvhoKIyGIiprkLuSuaqClGapZmRhpk+qSiiWGFT+pTmkuFS5ZZimLGo6AW4g4oorhAyC4IMsgmyjJyfn/wm/PMwKxw5pwzN+/36zWv13Du61zne65zz+H+3Mt1F6RSqVQAAADUcnVyXQAAAEAmCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDcAAEAiCDeQYAsWLIiCgoIoKCiICy+8cKf769q1axQUFETXrl3LXD5y5Mj06z3//PM1Xg/VV9lnxP/5/PPPY/To0XHIIYdEy5Yto06dOumxW716da7LyxvPP/98elxGjhxZZpvKvjsyobiGAQMG1Nhr1BYPPPBAejweeOCBXJcDWVUv1wVAbfDaa6/FX/7yl3j11VdjwYIFsXbt2mjQoEG0bt069t577+jbt28ceuihcdxxx0W7du22e/6CBQvS/8EMGDDAf74V+M1vfhOrV6+Oli1bxvDhw3NdDruozz77LI4++uh44403MtJfQUFBtZ/Tt2/fmDFjRkZen60b/AsWLIiIKDeEVdfzzz8fX/va16r9vGHDhsVvfvObjNQAlCbcQAXWrFkTF198cYwdO3a7ZUVFRfHZZ5/Fhx9+GJMnT47bb789CgoKYt26ddG0adNSbRcsWBA33nhj+t/CTfl+85vfxMKFC6NLly7CDTlzzz33pIPNfvvtF5dccknsvvvuUbdu3YiI7f7GyX8PPPBATJ48OSIyF26A/CPcQDk2b94cJ554YkyZMiUiIho0aBCnn356HHHEEdGxY8dIpVKxbNmymDFjRjz33HOxZMmSSKVSkUqlclx5zSne67mjunbtmujxqQ1Gjhxpw64KnnrqqYjYesRl4sSJsfvuu2es73HjxlWpXYsWLTL2mrm2s98dVZHr75Y+ffrETTfdVKW2e++9dw1XA7su4QbKceedd6aDTbdu3WLChAmxzz77lNk2lUrFa6+9FnfffXfUqeNSNqjtFi9eHBERHTp0yGiwiYgYNGhQRvsjP7Rt29ZnC3lAuIFyPPTQQ+nf77777nKDTcTWvbuHHXZYHHbYYdkoDahhGzdujIiIhg0b5rgSAKrDLmYox6xZs9K/H3300TvUR/EsQiUvOL3xxhvTs9iU/NnWokWL4ne/+12cffbZse+++0azZs2iQYMG0b59+xgwYEDcfPPNsWbNmmrXtGDBgrjiiiti3333jaZNm0br1q3jiCOOiLvvvju++OKLCp+7szMeVTRbWnHfCxcujIiIhQsXljlOxRMzHHLIIVFQUBANGjSI5cuXV/ray5cvj/r160dBQUEceuihO1T/tu9//fr1ceutt8Zhhx0W7du3jzp16pS553bz5s1x3333xde//vXYY489olGjRtGyZcvYf//94yc/+UmVT9lZvXp1/OxnP4v9998/mjVrFq1atYqDDjoofv3rX8f69evLrHFb1Zkt7Z133onLLrssevfuHS1atIjGjRvHXnvtFd/61rfi6aefrrTebWevWr9+fdxyyy1x0EEHRcuWLaNp06ax3377xYgRIzI++9hTTz0VQ4YMia5du0bjxo2jRYsW0adPn7j88svj3XffLfM5JWeYqmg9zNXsU5MnT466detGQUFB7LXXXhX+/c+fPz9atGgRBQUF0bx58/jggw9KLS9rPXjqqafi9NNPj86dO0fDhg2jc+fO8c1vfjNeffXVna69Ot8d//nPf+K//uu/4uCDD4527dpF/fr1o0WLFnHggQfGD3/4w/jXv/5V5ilo5c2WNmDAgCgoKEhfb1OybcmfXJ+uWVxnVSafqE7byuzs99OFF16YrqW47d/+9rc45ZRTonPnzlG/fv1o2bLlTtcJVZYCytSoUaNURKQiIrVw4cId6mPSpEnpPir72fZ5BQUFlT6nXbt2qRdffLHc158/f3667be//e3UxIkTUy1atCi3v4MPPjhVWFhYbn9dunRJRUSqS5cuZS6/4YYb0n1NmjSp0nrK6ruyn/vvvz+VSqVS9913X/qxm2++udyai40ePTrdfsyYMZW2L0vJ9z937txUz549t6vv9NNPL/WcN954I9WtW7cK31ODBg1S99xzT4WvPWPGjFSnTp3K7aNPnz6pRYsW7fRnlEqlUlu2bElde+21qTp16lRY96mnnppat25duTUXt+vfv3/qgw8+SPXq1avcvrp27ZpasGBBhWNQFWvXrk2dfPLJFdZdt27d1PXXX7/dc++///5qrYPVUd7fenVde+216X7OPffcMtts3rw5deihh1ZY77brwdChQ8t9v3Xq1EmNHDmy3JpKfs/dcMMNZbapbL0srnv48OGpunXrVvoZPP/889s9v+T6VlL//v2r9LmWV3tFSr73bV+3ukrWubNtS67LFa2vmfh++va3v51uO2vWrNRpp522XR8tWrSo9D1BpjgtDcrRvXv3eOeddyIi4o477ohf//rX1e5jv/32i3HjxsU777wT119/fUREDB48OM4999wKn7dhw4ZIpVLRp0+f+NrXvha9evWKNm3axIYNG2Lx4sXxxBNPxBtvvBErVqyIU089NWbMmFHpHtGFCxfG4MGDY+3atXH22WfHiSeeGE2aNIm33347xowZE4WFhTFt2rQYOHBgvPzyy+lZobLlD3/4Q3z22WdxySWXxIoVK6Jdu3bxhz/8Ybt2Bx54YEREnHvuuXHFFVfE6tWrY8yYMfFf//Vf5fadSqXi3nvvjYiI5s2bx+DBg3eq1o0bN8aZZ54Zs2bNiv79+8cZZ5wRnTp1iuXLl8e6devS7aZMmRLHHnts+qjKscceGyeffHLssccesWHDhnj11VfjwQcfjM8++ywuvfTSaNiwYZn3/1m2bFkcf/zxsWLFioiI6NmzZ1x44YXRtWvXKCwsjEcffTQmT54cgwcPjqKiop16bxERP/3pT9Pre7169WLIkCExYMCAaNiwYbz55ptx7733xqpVq2L8+PFxyimnxKRJkypcX9auXRsDBw6M2bNnx9e//vU4+eSTo3Xr1jFv3ry46667YvHixbFgwYK48MILY9KkSTtcd1FRUZx44onpIw2tW7eOiy66KA444IDYuHFjTJo0Kf7617/GF198EaNGjYqioqL45S9/mX7+Mccck77Yv6L1sHgdzIUbb7wx/vWvf8WUKVPib3/7W5x44onbrTM33nhjvPbaaxERcc4551R6T6nf/va38cQTT0Tbtm3j4osvjv333z8+++yzmDBhQjz22GOxZcuWGDlyZLRp0yYuu+yyGnlfqVQqzjrrrHjyyScjIqJu3boxaNCg+NrXvhbt27ePzz77LN5///2YOHFizJgxo1qTB9x0001RWFgY1113XfqoXVmTOvTs2TMzb6aWyNT3U0nDhw9PX596/vnnxz777BPr169PX78KWZHbbAX5a9SoUdvtkf/HP/6RWr16dbX7qsqezZIWLFiQevvttyts8/DDD6f3cF544YVltil5pCQiUvXq1Us98cQT27Vbvnx5qk+fPul2t956a5n91eSRm6q+RkmXX355pUchUqlU6t///ne63SWXXFJpv+XZ9ujS7bffXm7bdevWpds3bdo09dRTT5XZbs6cOak999wzFRGpJk2apD7++OPt2gwZMiT9mmeddVZq48aN27W56aabStW2o5/RSy+9lD5q2Lx589TLL7+8XZulS5eWOmo1evToMl+rZD0NGjRI/fOf/9yuzYoVK1Jdu3ZNt5s6dWqZfVXFL37xi3Q/vXv3Tn300UfbtXnhhRdSTZs2TR+ReOWVV8rsqzrrYVWUHIudNXfu3FTz5s1TEZFq1qxZas6cOellL774Yvp7Yc8990ytWrWqzD5KrgcRW4/8LV++fLt248aNS9WrVy+9fi5atGi7Npk4cnPzzTen+9hzzz0r/P57/fXXyzzKV/z88o6gVOfISFXV1iM3mfx+KnnkJmLrEcWyvqMgW4QbKMf69etTBx100HaH1wsKClL77LNP6pvf/Gbq9ttvT7377ruV9lXdcFNVxf+pNG7cOLVp06btlm8bbq6++upy+5o5c2Z6o2iPPfZIFRUVbdcm38LNu+++m+7vW9/6VrntSoaDadOmVdpveUqGm2984xsVtr3tttvSbR988MEK2/7rX/9Kt73ppptKLVu6dGl647JDhw6ptWvXltvPMcccs9Ph5vTTT08v//3vf1/ua82cOTNdV8eOHcvcmCm57v385z8vt68//OEP6XajRo0qt11FNm7cmGrfvn06xFf0d3n33XenX++MM84os01Nhpuq/lR0OtGDDz6YbnfwwQenNm3alFq1alV6Q7Ru3bqpF154odznl1wP6tWrl3rnnXfKbfvTn/403faaa67ZbvnOhpt169alWrdunQ7Ble3YKU+uw011fubPn79T9e1suMnU91MqVTrc7LnnnqnPPvus0vqhJplQAMrRpEmTeP755+Pyyy+PBg0apB9PpVIxZ86cePjhh+NHP/pR9OnTJ/r16xePPfZY1mssnp3t888/j7fffrvCtnXr1q3wppj77bdfnHjiiRGxdRrcadOmZazOmtK7d+846qijIiLisccei1WrVm3X5pNPPonHH388IiL69esXBx10UEZeu7LTcx588MGIiOjUqVN861vfqrDtMccck55u+Jlnnim17Kmnnkqfavad73wnmjdvXm4/l19+eaV1V2Tjxo3piQLatm0b3/nOd8ptu99++8XJJ58cEVtPm6voovO6detWOF7HHnts+vfyLvavzMsvvxwff/xxREQMHDgwevfuXW7b7373u9GmTZuI2Dq+mzZt2qHXzKXzzz8/vV5NmzYtrr/++vj+978fixYtioiIa6+9Nv23UZkTTzwx+vTpU+7y4cOHp087LP5byqSnn346Pvnkk4iIGDJkSHz5y1/O+GtQWqa+n7b13e9+Nxo3bpyZImEHueYGKtC0adO4/fbb4/rrr4/HHnssnnvuuXjttddiyZIlpdrNmDEjvvGNb8R3vvOdGDNmTMbudTNlypT4y1/+Eq+99lrMmzcv1q1bF5s3by6z7Ycffhhf+cpXyu2rT58+0aFDhwpf75hjjknfvHDatGk7PKtYNn3/+9+PF198MTZs2BB//vOf40c/+lGp5X/+859jw4YNERHxve99LyOvWa9evTjkkEPKXb527dp02OzUqVP6OoKKNGvWLCJKz9IXEfH666+nf992FqhtVba8MjNmzEhv6A8YMCDq169fYfvjjz8+/vnPf0bE1nW1f//+Zbbr0aNHtGrVqtx+vvSlL6V/LyugVsXUqVNL1VWRBg0aRP/+/ePxxx+PjRs3xowZM+KrX/3qDr3ujqjqTTwru7bnrrvuildeeSXmz58fN998c/rxQw89NH72s59VuZ6S4bIsHTt2jF69esU777wT//nPf2LNmjUZvcHoSy+9lP7961//esb6zbbq3MSzffv2NVxN+TL5/bStI488cucLhJ0k3EAVtGvXLi699NK49NJLIyLi448/jilTpsTTTz8dDz30UKxduzYiIu6///7Ye++9Y8SIETv1eps2bYqLL744/vznP1f5OcU1lKd79+6V9lGyzdKlS6v82rn0jW98I4YNGxYrV66MMWPGbBduiicSaNKkSQwZMiQjr9m6deto1KhRucsXLVoUW7ZsiYiI6dOnxxlnnFHlvov3YBf76KOP0r/vtddeFT63ZcuW0apVqx0OCCVfq6L7OhXr0aNH+vdly5aV265t27YV9lPyXjLFQbS6aqr2mpCpGz3utttu8de//jWOOuqo9NG95s2bx0MPPRT16lX9v/eqfje88847kUqlYtmyZRkNNx9++GH69169emWs32yrLTfxzOT307YyfcNb2BFOS4Md0L59+zjttNPirrvuirlz55baW3XzzTfH559/vlP9//CHP0wHm4YNG8YZZ5wRv/zlL+NPf/pTjB07NsaNGxfjxo0rdRpSZfeoadKkSaWv27Rp0/Tvn3766Q5Wn10lZ/CZOXNmqVl5XnvttZg5c2ZEbJ01KlMbZJWddrEj9x8qtu1sZyU/h5KfT3mq0qY8O/NaJWeJ21amjmRWpKZqz3df+tKXSr2Xgw46KLp161atPnL93VByx0zxEQJqTia/n7bllDTygXADO6lt27bx8MMPp/eUrlu3rtQpMtW1cOHC9NGGzp07x3vvvRePP/54XHPNNXHBBRfEWWedFYMGDYpBgwbFvvvuW+V+P/vss0rbFE8JGlG7NjK+//3vp29m98c//jH9eMnfM3VKWlWUHLsLL7wwUlsnb6nyT0klNyqr8hlWpU1V6q7u+lLRtUDZUJtr31FbtmyJ888/v9TG6qRJk+Luu++uVj+5/m7Ybbfd0r/Xlp0quVZ85GVHZPL7CfKRcAMZ0Llz51KnuezMKV0l77599dVXV7gXtvgu6lWx7R3KK2tTm04v2GeffeKYY46JiIi///3v8emnn8a6devikUceiYit58IffvjhWaun5NgV3ytpR3Xq1Cn9+/z58ytsu3r16kpPG6nqa82ZM6fS9iXblHxuLtTm2nfUr371q5g8eXJEbL1upvjI5E9+8pN4//33q9xPdb4bCgoKomPHjjtQbfk6d+6c/v29997LaN+1ScnTMyub5KKwsHCHXyeT30+Qj4QbyJCSM6ptu2ez5Gk5le35Wr58efr3vffeu8K2lc1cU9K7775bqu+ylLyB4sEHH1zlvjOpeKyqu4fw+9//fkRs3fP78MMPx8MPP5zeC5zNozYRW6/RKr524I033qhWCN1Wydndnn/++QrbFm/o7qgDDjggvR4///zzlZ6C8txzz6V/z+YF+WUp+fol6yrL5s2b02PVsGHD6Nu3b43WVhOmTJkSI0eOjIitG6t///vf46677oqIrbMnDhkypMqzwP373/+ucPmyZcvSYalHjx4Zvd4mIkrN6laVi9t3VHW+h3Oh5KQb205aU9Lq1avjP//5zw6/Tia/nyAfCTdQjsqCQEkLFy5MX9sREdtNq1oy7JQ8vaMsJc9/nzt3brnt/vGPf8Rbb71V5Rq/+OKLuP3228td/t5778XEiRMjYuue1FyFm+KxqmyctjVo0KD0HuU//vGP6VPSGjZsGOeff35mi6yCCy64ICK2bkRdffXVO9zPwIED06c83n///RWetnPHHXfs8OtEbB2rU045JSIiVqxYkZ4utizvvfde/O///m9EbJ1Nq3ha8lw54ogj0rMBjh8/PmbPnl1u2/vvvz+953vgwIGldkzUBuvWrYtvfetbUVRUFAUFBfGnP/0p2rRpE0OGDInzzjsvIrbOfHfNNddUqb8JEyZUeKTn9ttvT1/Td+aZZ+78G9jGySefHK1bt46IiIcffrjUd2kmVed7OBdKTl9eckfTtu66665Kr7GsTKa+nyAfCTdQjoMPPjguvvjiUlPxluXDDz+Ms846K/2fzWGHHbbdqWQlZ7maPn16hf2V3AN9yy23lDnz1dSpU+O73/1upe9hW7/+9a9j/Pjx2z2+YsWKOPfcc9N76kve1yLbisdq5cqV6Xt2VEX9+vXTYzJt2rT053bWWWelN5yy6bLLLos99tgjIiL+9re/xY9+9KPYuHFjue3XrFkTv/3tb7c76tCpU6c455xzImLrHvQLL7ywzD3yv/jFL+Jf//rXTtd95ZVXpq9f+vGPf1xqgoZiy5Yti7PPPju9vvz4xz/OeUBo0KBBDBs2LCK2Hpk5++yz0/e9KemVV16JK6+8MiK27sm/6qqrslpnJvzwhz9M7/j4yU9+Escdd1x62Z133pn+G/qf//mfePbZZyvtr6ioKAYPHhwrVqzYbtk///nPuOWWWyJi646XH/zgB5l4C6U0bdo0vYG9adOmOO200yoMODNmzNihow3V+R7OheL7RkVsPeWwrO/+p59+On7+85/v9Gtl6vsJ8pGpoKEcmzZtinvvvTfuvffe6N69exx99NFxwAEHRLt27aJOnTqxfPnyePXVV+OJJ55Iz47WrFmz9KkhJbVq1Sr69esXb775ZkyaNCkuvfTSOPbYY0tdyHzSSSdFxNZwdPDBB8e0adNiwYIF0bNnz7j00ktj3333jc8//zwmTZoUf/vb3yKVSsWQIUPir3/9a5Xez4ABA2LGjBnx9a9/Pc4+++w48cQTo0mTJvH222/HmDFj0hs2X/3qV9Mbiblw7LHHpk9NOfPMM+MHP/hBdOrUKX1KyZe//OVS90Up6Xvf+16MHj261MW2l1xySc0XXYZmzZrFE088EQMGDIh169bFHXfcEWPHjo1zzjkn+vbtG7vttlusW7cu5s+fH1OnTo1///vfsWnTpjKn/7711lvj2WefjRUrVsRjjz0Wffv2jQsvvDD22muvWLFiRTz66KMxefLkOOyww2LRokWxZMmSHZ6h7Igjjoif/OQnccstt8TatWvjyCOPjPPOOy8GDBgQDRo0iBkzZsSYMWPS1/YceeSR8ZOf/GSnxipTrrrqqvjnP/8Zr776asycOTN69+4dF198cfTt2zc2btwYzz//fDz00EPpUPbTn/40J/dyeuKJJ6rcduDAgaXuN/Twww+n15F+/frFL37xi1Ltd9ttt3jooYfiqKOOii+++CK+/e1vx9tvv13hdNyDBg2KJ554Ivr06RPf+9734stf/nJ89tlnMXHixHj00UfTp3DdfPPN6Q3iTLvyyivjpZdeiieffDIWLlwY/fr1izPOOCMGDBgQ7du3j88//zxmz54dzzzzTLz++usxadKk6NKlS7Ve49hjj00fvb7ooovixz/+cXTp0iW9I6d79+5Vmha7PIWFhVX+bFu0aBFf+9rXSj128MEHR//+/WPy5MnxwQcfxIEHHhiXXnppdOvWLT755JOYMGFC/OMf/4gePXpEkyZN4s0339zhWjP5/QR5JwWU6fjjj08VFBSkIqJKP3369ElNmzat3P6eeuqpVN26dct9fklz585NdenSpdy2DRs2TI0ZMyZ1//33px+7//77t3vN+fPnp5d/+9vfTk2cODHVokWLcvs9+OCDU4WFheW+h+KaunTpUubyG264Id3XpEmTKq2nLOvWrUv16NGj3BrLep8lnXzyyem2++yzT4Vtq6uy91+W9957L9W3b98qrUMNGzZMPf3002X28+abb6Y6depU4fq3cOHC1Je+9KVURKT233//Mvup7DNKpVKpLVu2pK655ppUnTp1Kqx34MCBqXXr1pX73ovb9e/fv9Jxqk7biqxZsyZ10kknVVh33bp1U9ddd12F/ezIZ12Rqn6PbPuzatWqdB/z589P//02adIk9f7775f7ejfeeGO6j9NPP3275duuB0OHDi23hoKCgtQNN9xQ7mtNmjQp3ba8dlUZz02bNqWGDh1a6XoXEanJkydv9/zK1qGioqLUkUceWW6fFb3Hqrz36vz07du3zP4WLFiQ6tq1a7nP23fffVP/+c9/Uv37908/VpbK/m8olonvp29/+9vpNvPnz6/2GEKmOS0NyvHMM8/EokWL4t57743vfOc78dWvfjXatWsXDRo0iPr160fr1q2jX79+8d3vfjeefPLJmDFjRqmLv7d18sknx8svvxxDhgyJvfbaq8L7AXTr1i2mT58e1157bfTu3TsaNWoUzZo1i3333Tcuu+yymD59elx00UXVfk8nnHBCzJgxI4YPH57e+9eiRYs47LDD4s4774xXXnkl2rRpU+1+M6lZs2bx2muvxYgRI+LAAw+MFi1aVOsoRMlTdLI9kUBZevXqFW+++WY88cQTccEFF8Q+++wTu+22W9StWzdatmwZffv2jQsuuCAeeOCB+Oijj9JH8LZ1wAEHxHvvvRfXXXdd9OnTJ/3ZHXjggXHzzTfHlClTYo899kifyrIzp+IVFBTEL3/5y5gxY0YMHTo0evbsGc2aNYtGjRpFly5d4txzz43//d//jfHjx+fdlOG77bZbPP300zF+/PgYPHhw7LnnntGwYcNo1qxZ9OzZM4YOHRpvvfVWjBo1KtelVssXX3wR5513Xnra59tuuy169uxZbvsRI0bEEUccERFbr8+75557Kuz/zjvvjPHjx8dpp50Wu+++ezRo0CB23333GDx4cLz88svpyQtqUv369ePOO++Mt956K4YNGxZf/vKXo2XLlum/la985Stx+eWXxwsvvBBHH310tfuvW7duPPvsszF69Og47LDDolWrVjk7/bY8Xbp0iTfffDNGjBgRvXv3jsaNG8duu+0W/fr1i1/96lfx+uuvV+kmtVWVqe8nyCcFqVQeThkCsIOOPPLIePnll6N+/frx4YcfRvv27XNdUta8++67sd9++0VExOWXX17hBBLs2kaOHBk33nhjRGy9eH3AgAG5LQggQxy5ARJj5syZ8fLLL0dExBlnnLFLBZuIKHW9l41VAHZFwg2QGD/72c/Sv//oRz/KYSWZ9+qrr8bmzZvLXX7vvfem70zfqVOnOO2007JVGgDkDbOlAbXWBx98EB988EGsXbs2HnvssfRMRccdd1z6eoOkuP766+Ott96KU045JQ466KDo2LFjFBUVxbx58+LJJ5+MqVOnptvec889pWbYAoBdRV6Fm/Hjx8esWbNi0aJFsWbNmti8eXO0bNkyevfuHaeffnqNTUEJ1E5/+ctf0tcNFGvdunWlF0/XVoWFhfHggw+We3PNRo0axT333BNf//rXs1wZAOSHvAo348aNiw0bNkSXLl1izz33jIiIxYsXxwsvvBCvvPJKXHXVVdGvX78cVwnkmzp16sSXvvSlOProo+PnP//5djdRTYLi2ayee+65mDdvXqxYsSLWrVsXLVu2jO7du8dxxx0XQ4cOjU6dOuW6VADImbyaLW3WrFnRrVu37e50/cwzz8SYMWOiVatWcffdd+/wzekAAIDkyquU0LNnz+2CTcTWe3N07NgxVq1aFUuXLs1BZQAAQL7Lq3BTkeKjNfXq5dWZdAAAQJ6oFeFm8uTJsXTp0ujUqdMud98KAACgavLyMMiTTz4Zixcvjo0bN8aSJUti8eLF0apVqxg2bJjrbQAAgDLlZbh56623YubMmel/t2nTJi6//PIqz4B0xRVXbPdY06ZNY9SoURmrEQAAyC95GW6uv/76iIhYv359LFq0KMaOHRsjR46Mc889N84888wcVwcAAOSjvJoKujxFRUVx3XXXxfz58+MXv/hFdO/efaf6W7ZsWbRp0yYitt4Uj5rXtm3biDDe2WTMs8t4Z5fxzj5jnl3GO/uMeXYVj3f9+vUz2m9eHrnZVr169eLwww+PefPmxRtvvLHT4aZknqsF2S5RjHf2GfPsMt7ZZbyzz5hnl/HOPmNeu9Waq/ObN28eERFr167NcSUAAEA+qjXh5r333ouIiA4dOuS4EgAAIB/lTbh5//3345VXXokvvvii1ONFRUXx9NNPxwsvvBANGjSIww8/PEcVAgAA+SxvrrlZvnx53HXXXdG8efPo1q1bNG/ePNatWxeLFi2KVatWRf369WPo0KHpi48AAABKyptw07t37zjjjDPivffei0WLFsXatWujXr160b59+zjkkEPilFNOiY4dO+a6TAAAIE/lTbhp3759fPOb38x1GQAAQC2VN9fcAAAA7AzhBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASIR6uS6g2MaNG+Ott96KN954I+bOnRsrVqyILVu2RMeOHeOQQw6JU089NRo1apTrMgEAgDyVN0duXnrppbjlllti0qRJkUqlom/fvtGzZ8/4+OOP45FHHolrrrkm1qxZk+syAQCAPJU3R27q1asXJ5xwQgwcODA6deqUfnzVqlUxevTomD9/fjzwwAMxbNiwHFYJAADkq7w5ctO/f/+4+OKLSwWbiIhWrVrFRRddFBERU6dOjaKiolyUBwAA5Lm8CTcV6dKlS0REbN68OdatW5fjagAAgHxUK8LN8uXLIyKibt260axZsxxXAwAA5KO8ueamIk899VRERBxwwAFRv379SttfccUV2z3WoEGDGD16dEREtG3bNurV2/rW27Vrl8FKKY/xzj5jnl3GO7uMd/YZ8+wy3tlnzLOreLwzLe+P3EyfPj0mTZoUdevWjcGDB+e6HAAAIE8VpFKpVK6LKM+HH34Y119/faxfvz4uvPDCOOWUUzLS70cffRRt27aNiIgVK1ZkpE8qVrwXxHhnjzHPLuOdXcY7+4x5dhnv7DPm2VU83lU5K6s68vbIzcqVK+OXv/xlrF+/Pk499dSMBRsAACCZ8jLcrF27Nm666aYoLCyMAQMGxPnnn5/rkgAAgDyXd+Hm888/j1/96lexZMmS+OpXvxqXXnppFBQU5LosAAAgz+VVuNm8eXP893//d8ydOzf69u0bw4cPjzp18qpEAAAgT+VNctiyZUv89re/jXfffTd69eoVV155ZY1NEQcAACRP3qSHCRMmxNSpUyMionnz5jFmzJgy251//vmx2267ZbM0AACgFsibcPPpp5+mfy8OOWU5++yzhRsAAGA7eRNuzjnnnDjnnHNyXQYAAFBL5c01NwAAADtDuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABJBuAEAABKhXq4LKGnevHnx9ttvxwcffBBz5syJVatWRf369eOhhx7KdWkAAECey6twM3bs2Hj99ddzXQYAAFAL5VW46dGjR3Tt2jX23nvv2HvvveOSSy7JdUkAAEAtkVfhZtCgQbkuAQAAqKVMKAAAACSCcAMAACSCcAMAACRCXl1zkylXXHHFdo81aNAgRo8eHRERbdu2jXr1tr71du3aZbW2XZXxzj5jXnP2vfv2arWf/YMf1VAluy7rd/YZ8+wy3tlnzLOreLwz3m+N9ApQS1U3uGSqTwEIAHZeQSqVSuW6iPKcc845NXITz48++ijatm0bERErVqzIaN+UrXgviPHOHmNeuSPHPpzrEir00je+mesS8pb1O/uMeXYZ7+wz5tlVPN7169fPaL+O3AC7jHwPM9vatl5hBwAqJtwAiVXbwkxlhB0AqJhwAyRK0gJNRUq+V0EHAIQbIAF2pUBTnuIxEHIA2JXlVbiZPn16PPbYY6UeKyoqihEjRqT/fdZZZ8WBBx6Y7dKAPCTUbM/RHAB2ZXkVbtauXRtz5swp9VgqlSr12Nq1a7NdFpBnhJqqcTQHgF1NXoWbAQMGxIABA3JdBpCHBJodJ+QAsKvIq3ADsC2hJnOEHACSrk6uCwAoj2BTM44c+7CxBSCRHLkB8o4N7+xwJAeApBFugLwh1OSGkANAUjgtDcg5p0nlB58BALWdcAPklA3q/CJoAlCbOS0NyAkb0PnNqWoA1EaO3ABZJ9jUHj4rAGoT4QbIKhvLtY/PDIDawmlpQFbYQK7dnKYGQG3gyA1Q4wSb5PBZApDPhBugRtkYTh6fKQD5SrgBaoQphZPNZwtAPhJugIyz4btrEGAByDfCDZBRNnZ3PT5zAPKFcANkjI3cXZfPHoB8INwAGWHjFusAALkm3AA7zUYtxawLAOSScAPsFBuzbMs6AUCuCDfADrMRS3msGwDkgnAD7BAbr1TGOgJAtgk3QLXZaKWqrCsAZJNwA1SLjVUAIF8JNwDUKIEYgGwRboAqs5HKjrLuAJANwg1QJTZO2VnWIQBqmnADVMpGKZliXQKgJgk3QIVsjAIAtYVwA0BWCcwA1BThBiiXjVBqinULgJog3ABlsvFJTbOOAZBpwg0AAJAIwg2wHXvUyRbrGgCZJNwApdjYJNuscwBkinADAAAkgnADpNmDTq5Y9wDIBOEGAABIBOEGiAh7zsk96yAAO0u4AWxUkjesiwDsDOEGAABIBOEGdnH2lJNvrJMA7CjhBgAASAThBnZh9pCTr6ybAOwI4QYAAEgE4QZ2UfaMAwBJI9wAkJcEcACqS7iBXZCNRgAgiYQbAPKWIA5AdQg3AABAIgg3sIuxJ5zaxjoLQFUJNwAAQCIINwAAQCIIN7ALcXoPtZV1F4CqEG4AAIBEEG4AAIBEEG5gF+G0Hmo76zAAlRFuAACARBBuAACARBBuYBfgdB6SwroMQEWEGwAAIBGEGwAAIBGEGwAAIBGEG0g41yiQNNZpAMoj3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3ECC7Xv37bkuAWqEdRuAsgg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3kFCmygUAdjXCDQC1kgAPwLbq5bqAbW3atCmeeOKJePnll6OwsDCaNWsWffv2jcGDB0ebNm1yXR4AAJCn8urIzaZNm2LUqFExduzY2LBhQxx00EHRpk2beP755+OnP/1pLFu2LNclAgAAeSqvjtyMGzcuZs+eHT169IjrrrsuGjVqFBER48ePjwcffDDuvvvuuPHGG3NcJQAAkI/y5shNUVFRTJgwISIiLrroonSwiYg49dRTo0uXLvH+++/HvHnzclUiAACQx/Im3MyaNSvWr18fHTp0iL322mu75YccckhERLz++uvZLg0AAKgF8ibcLFy4MCKizGATEdGtW7dS7QAAAErKm2tuCgsLIyLKnRGtdevWpdpV5IorrtjusQYNGsTo0aMjIqJt27ZRr97Wt96uXbsdqpfqMd5ATfCdkj2+x7PLeGefMc+u4vHOtLw5crNhw4aIiGjYsGGZy4uvwSluBwAAUFLeHLlJpVI7tbyk2267rcLlhYWF0bZt24iIWLFiRZX7ZccV7wUx3kAm+U7JHt/j2WW8s8+YZ1fxeNevXz+j/ebNkZvGjRtHRMTGjRvLXF78eMlZ1AAAAIrlTbgpPpKycuXKMpd/8sknpdoBAACUlDfhpkuXLhERMX/+/DKXF9/fprgdAABASXkTbnr27BlNmjSJ5cuXlxlwpkyZEhERBx54YLZLAwAAaoG8CTf16tWLk046KSIi7rvvvlKzoo0fPz4WLlwYPXv2jO7du+eqRAAAII/lzWxpERFnnnlmzJw5M2bPnh3Dhg2Lnj17RmFhYcyZMyeaN28eQ4cOzXWJAABAnsqbIzcRW2+0ecMNN8RZZ50VDRo0iGnTpsXHH38c/fv3j5tvvjk6duyY6xIByBOzf/CjXJcAQJ7JqyM3EVsDzuDBg2Pw4MG5LgUAAKhF8urIDZA59moDALsa4QYAAEgE4QYAAEgE4QYAAEgE4QYAAEgE4QYAAEgE4QaAWsdsgACURbiBBLMBCADsSoQbAAAgEYQbAAAgEYQbAAAgEYQbAGqVl77xzVyXAECeEm4g4WwIAgC7CuEGAABIBOEGAABIBOEGgFrDaZYAVES4gV2ADUIAYFcg3AAAAIkg3AAAAIkg3MAuwqlp1HbWYQAqI9wAAACJINwAkPcctQGgKoQb2IXYQAQAkky4AQAAEkG4ASCvOeIIQFUJN7CLsaEIACSVcAMAACSCcAO7IEdvqC2sqwBUh3ADAAAkgnADuyh7xMl31lEAqku4AQAAEkG4gV2YPePkK+smADtCuAEAABJBuIFdnD3k5BvrJAA7SrgBAAASQbgB7Cknb1gXAdgZwg0AAJAIwg0QEfaYk3vWQQB2lnADQM4JNgBkgnADpNnABABqM+EGKEXAIduscwBkinADAAAkgnADbMeedLLFugZAJgk3QJlsdFLTrGMAZJpwA0DWCTYA1AThBiiXDVAAoDYRboAKCThkmnUKgJoi3ACVsjFKpliXAKhJwg0AWSHYAFDThBugSmyYAgD5TrgBqkzAYUdZdwDIBuEGqBYbqVSXdQaAbBFugGqzsUpVWVcAyCbhBtghNlqpjHUEgGwTboAdZuOV8lg3AMgF4QbYKTZi2ZZ1AoBcEW6AnWZjlmLWBQBySbgBMsJGLdYBAHJNuAEyxsbtrstnD0A+EG6AjLKRu+vxmQOQL4QbIONs7O46fNYA5BPhBqgRNnqTz2cMQL4RboAaY+M3mV76xjd9tgDkJeEGqFE2hJPFZwlAPhNugKywUVz7+QwByHfCDZA1No5rL58dALVBvVwXAOxaijeSjxz7cI4roSqEGgBqE0dugJyw0Zz/fEYA1DaO3AA54yhOfhJqAKitHLkBcs7GdP7wWQBQmzlyA+QFR3FyS6gBIAmEGyCvCDnZJdQAkCTCDZCXhJyaJdQAkETCDZDXhJzMEmoASDLhBqgVhJydI9QAsCvIi3CzYcOGmDp1anzwwQcxZ86cWLhwYRQVFcWQIUNi0KBBuS4PyCNCTvUINQDsSvIi3Cxbtix+97vf5boMoBYRciom1ACwK8qLcNOoUaM45phjonv37rH33nvHlClT4vHHH891WUAtUHIjflcPOgINALu6vAg3HTt2jEsvvTT972nTpuWwGqC22hWDjkADAP8nL8INQKYlOegINABQNuEGSLxtw0BtCzvCDABUjXAD7HLyOewIMgCw4xIZbq644ortHmvQoEGMHj06IiLatm0b9eptfevt2rXLam27KuOdfca86mb/4EflLtv37tuz/ppUzvqdfcY8u4x39hnz7Coe74z3m4lObr311li8eHG1nnPZZZdF9+7dM/HyADWmqiGk+Eu6qKioJssBACqQkXCzYsWKWLp0abWes3Hjxky8dJluu+22CpcXFhZG27ZtI2Jr7dS84r0gxjt7jHl2Ge/sMt7ZZ8yzy3hnnzHPruLxrl+/fkb7zUi4KT7dCwAAIFfq5LoAAACATBBuAACARBBuAACARMibqaB//etfx+rVqyMiYuXKlRERMXHixJg2bVpERLRs2TKuuuqqXJUHAADkubwJNwsWLNhudoqVK1emg445xwEAgIrkTbi58847c10CAABQi7nmBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASIR6uS4gImLJkiUxbdq0ePvtt+Ojjz6KNWvWRNOmTWPfffeNgQMHRq9evXJdIgAAkOfyItyMGjUqPvnkk2jcuHHss88+0aNHj/jwww9j6tSpMW3atLjgggti4MCBuS4TAADIY3kRbjp37hznnXdeHHrooVGv3v+V9Oyzz8Yf//jH+POf/xx9+/aNzp0757BKAAAgn+XFNTfXXXddHHnkkaWCTUTE8ccfH3379o0tW7bEq6++mqPqAACA2iAvwk1FunTpEhERq1atynElAABAPsv7cLN8+fKIiGjZsmVuCwEAAPJaXoebZcuWxfTp0yMi4qCDDspxNQAAQD4rSKVSqVwXUZYvvvgibrzxxpg1a1YcfvjhMXz48Co/94orrtjusQYNGsTo0aMjImLz5s3p63uKiooyUi8VM97ZZ8yzy3hnl/HOPmOeXcY7+4x5dhWPd0FBQWb7zUQnt956ayxevLhaz7nsssuie/fu5S6/7777YtasWdGhQ4e4+OKLd7ZEAAAg4TISblasWBFLly6t1nM2btxY7rKxY8fGs88+Gy1atIgRI0ZEs2bNqtX3bbfdVuHywsLCaNu2bURsrZ2a165du4gw3tlkzLPLeGeX8c4+Y55dxjv7jHl2FY93/fr1M9pvRsJN8elemTBhwoR45JFHokmTJjFixIjo2LFjxvoGAACSK68mFHjxxRfj/vvvj4YNG8bVV18dXbt2zXVJAABALZE34Wb69Olx1113Rd26dePKK6+Mnj175rokAACgFsmLcDNr1qz0dTLDhw+Pvn375rgiAACgtsnINTc76+abb45NmzZF+/btY9q0aTFt2rTt2vTs2TOOPfbYHFQHAADUBnkRbtavXx8RER9//HF8/PHH5bYTbgAAgPLkRbh55JFHcl0CAABQy+XFNTcAAAA7S7gBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASQbgBAAASoV6uC4iIWLhwYTz33HMxb968KCwsjHXr1kX9+vWjc+fOcdRRR8Xxxx8fdevWzXWZAABAHsuLcPP+++/HxIkTo127dtG5c+do3rx5rFu3LmbNmhVz5syJqVOnxrXXXhv16uVFuQAAQB7Ki7TQr1+/uOOOO6JDhw6lHl+9enWMGjUq3nnnnfj3v/8dJ5xwQo4qBAAA8l1eXHPToUOH7YJNRETLli1j0KBBERHxzjvvZLkqAACgNsmLcFOROnW2luiUNAAAoCJ5nRg+/fTTGD9+fERsPXUtUwoKCsr8nZpnvLPPmGeX8c4u4519xjy7jHf2GfParSCVSqVyXUSxjz76KB5//PFIpVKxZs2amD17dmzYsCGOO+64+N73vlflle2KK67Y7rGmTZvGqFGjMl0yAACQJ/LqyM2aNWti8uTJpR476aST4txzz5WiAQCACmXkyM2tt94aixcvrtZzLrvssujevXuZy7Zs2RKFhYUxderUePTRR6Nly5YxYsSIaN++/c6Wmnb11VdHRMTo0aMz1iflM97ZZ8yzy3hnl/HOPmOeXcY7+4x5dtXUeGfkyM2KFSti6dKl1XrOxo0by11Wp06daN++fZx66qnRvn37uOWWW+K+++5LD0ImbNq0KWN9UTnjnX3GPLuMd3YZ7+wz5tllvLPPmGdXTY13RsJNTSbcgw8+OBo1ahQzZsyIoqIis6YBAABlyvupoAsKCqJZs2axZcuW+PTTT3NdDgAAkKfyPtwsX748Vq5cGY0bN47ddtst1+UAAAB5Ki/CzT/+8Y9Yvnz5do8vXbo0br/99kilUtG/f//0DT0BAAC2lRf3ufnhD38YhYWF0bVr1+jQoUNEbJ2kYN68eZFKpaJXr15xzTXXRKNGjXJcKQAAkK/yIty8+OKL8eabb8bcuXNj9erVsWnTpmjWrFl07do1jjjiiDj66KMdtQEAACqUF+EGAABgZzkcAgAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJEK9XBeQz8aOHRuPPPJIREQMGzYsjjjiiBxXlBwLFy6M5557LubNmxeFhYWxbt26qF+/fnTu3DmOOuqoOP7446Nu3bq5LjMxlixZEtOmTYu33347Pvroo1izZk00bdo09t133xg4cGD06tUr1yUmzoYNG2Lq1KnxwQcfxJw5c2LhwoVRVFQUQ4YMiUGDBuW6vFpr06ZN8cQTT8TLL78chYWF0axZs+jbt28MHjw42rRpk+vyEmXevHnx9ttvp9fhVatWRf369eOhhx7KdWmJtHHjxnjrrbfijTfeiLlz58aKFStiy5Yt0bFjxzjkkEPi1FNPdb+/GjB+/PiYNWtWLFq0KNasWRObN2+Oli1bRu/eveP000+PPfbYI9clJtann34aw4cPj7Vr18buu+8ev/nNbzLSr3BTjqVLl8a4ceOioKAgzJadee+//35MnDgx2rVrF507d47mzZvHunXrYtasWTFnzpyYOnVqXHvttVGvnlU0E0aNGhWffPJJNG7cOPbZZ5/o0aNHfPjhhzF16tSYNm1aXHDBBTFw4MBcl5koy5Yti9/97ne5LiNRNm3aFKNGjYrZs2dHq1at4qCDDooVK1bE888/H9OnT4+bbropOnbsmOsyE2Ps2LHx+uuv57qMXcZLL70Uv//97yMiYo899oi+ffvG559/Hv/5z3/ikUceiZdffjlGjhwZLVq0yHGlyTJu3LjYsGFDdOnSJfbcc8+IiFi8eHG88MIL8corr8RVV10V/fr1y3GVyfSnP/0p1q1bl/F+bTmWIZVKxe9///to0qRJ7LPPPr7ca0C/fv3ijjvuiA4dOpR6fPXq1TFq1Kh455134t///neccMIJOaowWTp37hznnXdeHHrooaUC47PPPht//OMf489//nP07ds3OnfunMMqk6VRo0ZxzDHHRPfu3WPvvfeOKVOmxOOPP57rsmq1cePGxezZs6NHjx5x3XXXpfdijx8/Ph588MG4++6748Ybb8xxlcnRo0eP6Nq1a+y9996x9957xyWXXJLrkhKtXr16ccIJJ8TAgQOjU6dO6cdXrVoVo0ePjvnz58cDDzwQw4YNy2GVyXPVVVdFt27dokGDBqUef+aZZ2LMmDFxzz33xN133+1m8hk2c+bMmDx5chx33HHx3HPPZbRvn1QZ/vWvf8X7778fF1xwQTRt2jTX5SRShw4dtgs2EREtW7ZMn7LzzjvvZLmq5LruuuviyCOP3O5I2PHHHx99+/aNLVu2xKuvvpqj6pKpY8eOcemll8Zxxx0Xe+21l/8Yd1JRUVFMmDAhIiIuuuiiUqfnnHrqqdGlS5d4//33Y968ebkqMXEGDRoU55xzTnzlK1+Jli1b5rqcxOvfv39cfPHFpYJNRESrVq3ioosuioiIqVOnRlFRUS7KS6yePXtuF2wiIk444YTo2LFjrFq1KpYuXZqDypJr06ZN8cc//jE6d+4cp512Wsb797/tNlavXh0PPfRQfPnLX46jjjoq1+Xskoo3Ap2Slh1dunSJiK17ByFfzZo1K9avXx8dOnSIvfbaa7vlhxxySESEI+0kUvH39ObNm2vkNB7KZnukZjz66KOxfPnyuPjii2vk+mrhZhv33XdfbNq0KS6++OJcl7JL+vTTT2P8+PEREc5xzZLly5dHRNgzS15buHBhRESZwSYiolu3bqXaQZIUf0/XrVs3mjVrluNqdg2TJ0+OpUuXRqdOnaJ9+/a5LicxFi5cGOPHj48BAwZE7969a+Q1RNES3njjjXjttdfinHPO2e6wMDXjo48+iscffzxSqVSsWbMmZs+eHRs2bIjjjjsujjzyyFyXl3jLli2L6dOnR0TEQQcdlONqoHyFhYUREeXOiNa6detS7SBJnnrqqYiIOOCAA6J+/fo5riaZnnzyyVi8eHFs3LgxlixZEosXL45WrVrFsGHDnFacIVu2bElf037eeefV2OsIN//fhg0bYsyYMdGpU6c4/fTTc13OLmPNmjUxefLkUo+ddNJJce6550ZBQUGOqto1fPHFF3HXXXfF5s2b4/DDD0/v+YZ8tGHDhoiIaNiwYZnLi6/BKW4HSTF9+vSYNGlS1K1bNwYPHpzrchLrrbfeipkzZ6b/3aZNm7j88sv935hBEyZMiA8++CCGDh0azZs3r7HXSUy4ufXWW2Px4sXVes5ll10W3bt3j4iIv/71r7Fy5cr42c9+Zq9IFezseBfr2bNnPPLII7Fly5YoLCyMqVOnxqOPPhpvvfVWjBgxwqHg/y9T413SfffdF7NmzYoOHTo4DbMMNTHm7LjKpuQ3ZT9J9OGHH8Ydd9wRqVQqzj///OjatWuuS0qs66+/PiIi1q9fH4sWLYqxY8fGyJEj49xzz40zzzwzx9XVfoWFhfG3v/0tevfuHQMGDKjR10pMuFmxYkW1Z7PYuHFjRER88MEHMXHixDj66KNjv/32q4nyEmdnxrssderUifbt28epp54a7du3j1tuuSXuu+++uPrqq3e21ETI9HiPHTs2nn322WjRokWMGDHCOdxlyPSYs3MaN24cEeWPcfHjbnJIUqxcuTJ++ctfxvr16+PUU0+NU045Jdcl7RKaNm0avXr1imuuuSauu+66+Pvf/x7777+/HVc7acyYMVFUVJSVnamJCTejR4/e4edOnz49UqlULFq0KEaOHFlq2ZIlSyLi/zYGDz300DjppJN2ptRE2JnxrszBBx8cjRo1ihkzZkRRUZFZSiKz4z1hwoR45JFHokmTJjFixAg3PSxHTa7jVF/btm0jYusGX1k++eSTUu2gNlu7dm3cdNNNUVhYGAMGDIjzzz8/1yXtcurVqxeHH354zJs3L9544w3hZidNnz49mjZtGmPGjCn1+ObNmyNi65Gd4m3wq6++eqd2VNlqLGHBggXlLluyZEksWbLEIeEsKCgoiGbNmkVhYWF8+umnZvHKoBdffDHuv//+aNiwYVx99dXWZ2qN4qlw58+fX+by4vvbFLeD2urzzz+PX/3qV7FkyZL46le/GpdeeqlrUHOk+LqQtWvX5riSZFi/fn289957ZS7btGlTetkXX3yxU68j3ETEOeecE+ecc06Zy+68886YPHlyDBs2LI444ogsV7ZrWr58eaxcuTIaN24cu+22W67LSYzp06fHXXfdFXXr1o0rr7wyevbsmeuSoMp69uwZTZo0ieXLl8f8+fO3mxJ6ypQpERFx4IEH5qI8yIjNmzfHf//3f8fcuXOjb9++MXz4cDN15VDxxnZZNx2neh555JEyH//444/jsssui9133z1+85vfZOS1/MWQE//4xz/S8/aXtHTp0rj99tsjlUpF//79falnyKxZs+K2226LiIjhw4dH3759c1wRVE+9evXSpwTfd999pWZFGz9+fCxcuDB69uzp1BFqrS1btsRvf/vbePfdd6NXr15x5ZVXOi27hr3//vvxyiuvbHekoKioKJ5++ul44YUXokGDBnH44YfnqEJ2hL8acuKZZ56Jv/71r9G1a9f0HpEVK1bEvHnzIpVKRa9evWLIkCE5rjI5br755ti0aVO0b98+pk2bFtOmTduuTc+ePePYY4/NQXXJ9etf/zpWr14dEf93rcjEiRPT49+yZcu46qqrclVerXPmmWfGzJkzY/bs2TFs2LDo2bNnFBYWxpw5c6J58+YxdOjQXJeYKNOnT4/HHnus1GNFRUUxYsSI9L/POussR8syZMKECTF16tSI2Ho61LbXJhQ7//zzndWQIcuXL4+77rormjdvHt26dYvmzZvHunXrYtGiRbFq1aqoX79+DB061LV8tYxwQ06ce+658eabb8bcuXPjrbfeik2bNkWzZs1i//33jyOOOCKOPvpoR20yaP369RGx9fDvxx9/XG474SazFixYECtWrCj12MqVK9NBp127drkoq9Zq0KBB3HDDDTFu3Lh46aWXYtq0adG0adPo379/DB482AZIhq1duzbmzJlT6rFUKlXqMdciZM6nn36a/r045JTl7LPPFm4ypHfv3nHGGWfEe++9F4sWLYq1a9dGvXr1on379nHIIYfEKaecYtKdWqgg5eYAAABAAtg1DgAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJIJwAwAAJML/A7kjsyWt655fAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Expression\n",
    "expr = sy.simplify((y + h*ff(t, y))/y)\n",
    "\n",
    "print(expr)\n",
    "\n",
    "# Sig\n",
    "z = sy.Symbol('z')\n",
    "sigf = sy.lambdify(z, expr.subs(h*lam, z), 'numpy')\n",
    "\n",
    "# Make grid\n",
    "xx = np.linspace(-3, 3, 201)\n",
    "yy = np.linspace(-3, 3, 201)\n",
    "X, Y = np.meshgrid(xx, yy)\n",
    "z = X + Y*1j\n",
    "\n",
    "# Amplication factor of Explicit Euler (z=lambda h)\n",
    "sig = sigf(z)\n",
    "\n",
    "# Same scale of x and y\n",
    "plt.axis('equal')\n",
    "\n",
    "# Stability region\n",
    "plt.contourf(X,Y,abs(sig), levels=[0, 1])\n",
    "\n",
    "# Title\n",
    "plt.title('Stability region of Explicit Euler')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0130aacf",
   "metadata": {},
   "source": [
    "### 2nd-order Runge-Kutta Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "04b97b69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5*h*lambda*(h*lambda + 2) + 1\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Stability region of RK2')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzcAAAKOCAYAAACbY+fTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAABH6UlEQVR4nO3deXhV5bk34CcQAjLIPLUoKMikFfVonYUijlBnRWm1tloP9TjV6lct2mqxFk/VtlqlAw5HP2uPxaGWOtQBcRYqIjjAh8yCYIJAAGWI7O8PTvZJIAkJ7Oy9s7jv6+K6kr3WfteTd6+E9VvvWu8qSKVSqQAAAGjgGuW6AAAAgEwQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQboCdxvz586OgoCAKCgri/PPP3+H2evToEQUFBdGjR48ql99www3p7b300kv1Xg91t63PiP/1xRdfxJgxY+Lggw+ONm3aRKNGjdJ9t3LlylyXBxAREYW5LgBoeN588834v//3/8Ybb7wR8+fPj9LS0igqKop27dpFz549Y8CAAXHIIYfEkCFDomPHjlu9f/78+XH//fdHRMSgQYNi0KBB2f0BGpDf/OY3sXLlymjTpk1cccUVuS6HndTnn38eRx11VLz99tsZaa+goKDaZS1btowOHTrEvvvuG8OGDYsRI0ZEixYtamxv/vz5sccee0RERPfu3WP+/Pk1rr9y5co48cQT44033oiIiD59+sTzzz8f3bp1S68za9aseOaZZ+KVV16J6dOnx5IlS2Ljxo3Rtm3b2GeffeK4446L733ve9G+ffta/tRANgg3QK2tWrUqLrzwwhg/fvxWy8rKyuLzzz+Pjz/+OCZNmhR33HFHFBQUxOrVq7c6MJk/f37ceOON6e+Fm+r95je/iQULFkT37t2FG3Lm97//fTrY7LPPPnHRRRfFV77ylWjcuHFExDbDR12sWbMm1qxZE/Pnz48nn3wybrrppvjLX/4Shx56aEbaLy4ujmOPPTamTZsWERH77bdfPPvss9GpU6f0Ovvvv396+ZaWLVsWy5YtixdeeCFuvvnmGDt2bJx99tkZqQ3YccINUCsbN26M4447Lt56662IiCgqKoqTTz45Dj/88OjSpUukUqlYunRpTJs2LZ5//vlYvHhxpFKpSKVSOa68/mzr7PC29OjRI9H90xDccMMNccMNN+S6jLz31FNPRcTmEZdnn302vvKVr2Ss7ccff7zS96WlpfHOO+/Egw8+GMuXL4+FCxfGiSeeGNOmTYvu3bvv0LYWL14cQ4YMiZkzZ0ZExKGHHhpPPfVUtGnTptJ6M2bMiIiIRo0axeGHHx4DBw6Mnj17RvPmzWPevHnx5z//OaZPnx4rV66MESNGxKZNm2LEiBE7VBuQGcINUCt33XVXOtjsueee8cwzz8Ree+1V5bqpVCrefPPNGDt2bDRq5NY+aOgWLVoUERGdO3fOaLCJiDjllFO2eu28886La6+9NgYOHBgzZ86MlStXxk033RR/+tOftns7c+fOjaOPPjp9UuLoo4+Ov/3tb1WOOu26664xcuTIGDlyZOy+++5bLb/66qvj6quvjttvvz1SqVRccsklccIJJ0Tbtm23uz4gMxx1ALXy0EMPpb8eO3ZstcEmYvPZ3UMPPTQeeOCBaN68eTbKA+rR+vXrIyKiadOmWdtmp06d4tZbb01///e//3272/rggw/iyCOPTAebk046Kf7xj39UezndnDlz4uabb64y2ERsHtG59dZb46CDDoqIiBUrVsSTTz653fUBmSPcALVSfhlHRMRRRx21XW289NJLUVBQEN/4xjfSr914443pGZcq/tvSwoUL43e/+12ceeaZ0adPn2jZsmUUFRVFp06dYtCgQXHLLbfEqlWr6lzT/Pnz48orr4w+ffpEixYtol27dnH44YfH2LFj48svv6zxvduaLa02265utrTythcsWBAREQsWLKiyn8onZjj44IOjoKAgioqKYtmyZdvc9rJly6JJkyZRUFAQhxxyyHbVv+XPv3bt2rjtttvi0EMPjU6dOkWjRo2qPCu/cePGuPfee+Okk06K3XbbLZo1axZt2rSJfffdN370ox/V+nK/lStXxk9/+tPYd999o2XLltG2bds48MAD41e/+lWsXbu2yhq3VJfZ0t5777245JJLon///tG6devYZZddYo899ohvfetb8fTTT2+z3vLtlN9jtnbt2rj11lvjwAMPjDZt2kSLFi1in332iVGjRmV89rGnnnoqRowYET169IhddtklWrduHXvvvXdceuml8f7771f5nvvvvz9dc037Yfk+WB+OPPLI9NfLli2L0tLSOrcxderUGDhwYCxZsiQiIs4555x49NFHawxqtRmBKSgoiNNPPz39/fTp0+tcG1APUgC10KxZs1REpCIitWDBgu1qY+LEiek2tvVvy/cVFBRs8z0dO3ZMvfLKK9Vuf968eel1v/Od76SeffbZVOvWratt76CDDkqVlJRU21737t1TEZHq3r17lct/9rOfpduaOHHiNuupqu1t/bvvvvtSqVQqde+996Zfu+WWW6qtudyYMWPS648bN26b61el4s8/Z86cVN++fbeq7+STT670nrfffju155571vgzFRUVpX7/+9/XuO1p06alunbtWm0be++9d2rhwoU7/BmlUqnUpk2bUj/5yU9SjRo1qrHuYcOGpVavXl1tzeXrDRw4MPXRRx+l+vXrV21bPXr0SM2fP7/GPqiN0tLS1AknnFBj3Y0bN05df/31W733vvvuq9M+WBfV/a5vad26dZXWXbp0aZXrVfxdqvhZv/rqq5V+xy+66KLUl19+Wed6q3PXXXel2x45cmTG2gW2n3tugFrp1atXvPfeexERceedd8avfvWrOrexzz77xOOPPx7vvfdeXH/99RERMXz48G3ONLRu3bpIpVKx9957xze+8Y3o169ftG/fPtatWxeLFi2KJ554It5+++0oLi6OYcOGxbRp07Y5mrJgwYIYPnx4lJaWxplnnhnHHXdcNG/ePKZPnx7jxo2LkpKSmDJlSgwdOjRee+219KxQ2fLHP/4xPv/887jooouiuLg4OnbsGH/84x+3Wu+AAw6IiIizzz47rrzyyli5cmWMGzcu/s//+T/Vtp1KpeKee+6JiIhWrVrF8OHDd6jW9evXx2mnnRYzZ86MgQMHxqmnnhpdu3aNZcuWxerVq9PrvfXWW3H00UenR1WOPvroOOGEE2K33XaLdevWxRtvvBEPPPBAfP755zFy5Mho2rRplc//Wbp0aRxzzDFRXFwcERF9+/aN888/P3r06BElJSXx17/+NSZNmhTDhw+PsrKyHfrZIiJ+/OMfp/f3wsLCGDFiRAwaNCiaNm0a77zzTtxzzz2xYsWKmDBhQpx44okxceLEGveX0tLSGDp0aMyaNStOOumkOOGEE6Jdu3Yxd+7cuPvuu2PRokUxf/78OP/882PixInbXXdZWVkcd9xx6emO27VrFxdccEHst99+sX79+pg4cWL8+c9/ji+//DJGjx4dZWVlcfPNN6ffP3jw4PTN/jXth+X7YH2oOKrUtGnTSjOabctzzz0Xp5xySnz++ecREfGjH/2o0mVumVD+NzEidniyAyBDcp2ugIZh9OjRW52R/9vf/pZauXJlnduqOILzs5/9bJvrz58/PzV9+vQa13n44YdTjRs3TkVE6vzzz69ynYpndyMiVVhYmHriiSe2Wm/ZsmWpvffeO73ebbfdVmV79TlyU9ttVHTppZducxQilUqlXnzxxUpnsrfXlqNLd9xxR7Xrrl69Or1+ixYtUk899VSV682ePTu1++67pyIi1bx589Snn3661TojRoxIb/P0009PrV+/fqt1brrppkq1be9n9Oqrr6ZHDVu1apV67bXXtlpnyZIllUatxowZU+W2KtZTVFSU+vvf/77VOsXFxakePXqk15s8eXKVbdXGL37xi3Q7/fv3T33yySdbrfPyyy+nWrRokYqIVKNGjVKvv/56lW3VZT+sjYp9UZMzzjij0ohXdbYcuXniiSdSTZs2Tb92ww03ZKTuilasWJFq165dehvvvvtuxrcB1J1wA9TK2rVrUwceeOBWl6QUFBSk9tprr9Q555yTuuOOO1Lvv//+Ntuqa7ipre985zupiEjtsssuqQ0bNmy1fMtwc80111Tb1owZM9JhabfddkuVlZVttU6+hZv3338/3d63vvWtaterGA6mTJmyzXarUzHcnHHGGTWue/vtt6fXfeCBB2pc94UXXkive9NNN1VatmTJklRhYWEqIlKdO3dOlZaWVtvO4MGDdzjcnHzyyenlf/jDH6rd1owZM9J1denSpcrAVXHf+/nPf15tW3/84x/T640ePbra9Wqyfv36VKdOndIhvqbfy7Fjx6a3d+qpp1a5TjbDTWlpaWrSpEmpb37zm5XWe/rpp6ttr+LvUtOmTdOfRUFBQerXv/51Rmre0ve///30NocOHVov2wDqzoQCQK00b948Xnrppbj00kujqKgo/XoqlYrZs2fHww8/HJdddlnsvffesf/++8ejjz6a9RrLH/L3xRdfbPPm3saNG9f4UMzyJ5BHbJ4Gd8qUKRmrs770798/fQP2o48+GitWrNhqnc8++ywee+yxiNj8oMIDDzwwI9u+5JJLalz+wAMPRERE165d41vf+laN6w4ePDg93fA///nPSsueeuqp9KVm3/3ud6NVq1bVtnPppZdus+6arF+/Pj1RQIcOHeK73/1utevus88+ccIJJ0TE5svmyi8Fq0rjxo1r7K+jjz46/XV1N/tvy2uvvRaffvppREQMHTo0+vfvX+263/ve96J9+/YRsbl/N2zYsF3b3F5bTlCw6667xsCBAyvNjnb77bfH8ccfX6v21q9fn95Hdt999zjvvPMyXvODDz6YnpZ61113jTvuuCPj2wC2j3AD1FqLFi3ijjvuiI8//jjGjh0bp59+enz1q1/dar1p06bFGWecEd/73vdi06ZNGdv+W2+9FZdeemkcdNBB0b59+ygqKqp0UDRy5Mj0uh9//HGNbe29997RuXPnGtcZPHhw+uuGEG4iIv793/89Ijbfp/Tggw9utfzBBx+MdevWRUTE97///Yxss7CwMA4++OBql5eWlqbDZteuXePJJ5+MJ554osZ/LVu2jIjKs/RFRPzrX/9Kf10+61h1trV8W6ZNm5Y+0B80aFA0adKkxvWPOeaY9Nflz4SqSu/evWucjavi71RVAbU2Jk+eXGVdVSkqKoqBAwdGxOZgMG3atO3aZn3Yf//9Y8aMGfHDH/6w1u9p06ZNug8XLFgQQ4YM2e5+rMorr7wSF110UURsDmb33HNP7LnnnhlrH9gxJhQA6qxjx47pB9xFRHz66afx1ltvxdNPPx0PPfRQerrW++67L3r27BmjRo3aoe1t2LAhLrzwwioP1quzrSlje/Xqtc02Kq5TPo1svjvjjDPi8ssvj+XLl8e4cePisssuq7S8fCKB5s2bZ+yJ6u3atYtmzZpVu3zhwoXpkDt16tQ49dRTa932Z599Vun7Tz75JP31HnvsUeN727RpE23btt3uA9uK26rpuU7levfunf566dKl1a7XoUOHGtupOEVxeRCtq/qqvT6UT1oQsXnUdf78+fHQQw/F+++/H++880787ne/i7vvvrvWDwRu3bp1PPfcczFo0KBYsmRJvPPOO3HMMcfE888/H23atNmhWv/1r3/FsGHD0p/Lr3/96zjjjDN2qE0gs4zcADusU6dO8c1vfjPuvvvumDNnThxxxBHpZbfcckt88cUXO9T+f/zHf6SDTdOmTePUU0+Nm2++Of7rv/4rxo8fH48//ng8/vjjlS5D2tYzamrzcNGKD/hbs2bNdlafXRVnGJsxY0alEYQ333wzZsyYERERZ511VrRu3Toj29xll11qXL49zx8qt+VsZxU/h+oewFhRbdapzo5sq+IscVuq7UH6jqiv2uvDKaeckv53zjnnxLXXXhszZsxI/z7/4Q9/SM+uWFt77bVXTJw4Mbp27RoREW+//XYce+yxO7QvTp8+PY477rj0iZObbropLr/88u1uD6gfwg2QUR06dIiHH344Cgs3DwyvXr260iUydbVgwYL0aEO3bt3igw8+iMceeyyuvfbaOO+88+L0009PHxj16dOn1u2WTw9bk/IpiyMifZlUQ/Dv//7v6Qehlt8XsOXXmbokrTYq9t35558fqc2T2dT6X0UVD8Jr8xnWZp3a1F3X/aWme4GyoSHXHrH5cq9f//rX6XvCxowZU+e/I717946JEydGly5dImLzpaUVw0ldvP/++zFkyJD0SOL111+/wyPSQP0QboCM69atW6XLXHbkkq4XXnghfYB7zTXX1Hhte/lT1Gvjo48+qtM65Te4NwR77bVX+n6h//7v/441a9bE6tWr45FHHomIzfcbHXbYYVmrp2LfVXwuyPYoPxMfETFv3rwa1125cuVWl7Vt77Zmz569zfUrrlPxvbnQkGsv17hx47j99tsjImLTpk3xox/9qM5t9OnTJ1588cX0/XVvvfVWHHfccXUanZo5c2YcffTR6ecqXXPNNfHzn/+8zrUA2SHcAPWi4oxqW456VLwsZ8sz81tatmxZ+uuePXvWuO6WM2vV5P3336/UdlUqPkDxoIMOqnXbmVTeV9vqpy2VTyywZs2aePjhh+Phhx9OX6qUzVGbiM33aPXr1y8iNl8eVJcQuqWKs7u99NJLNa47adKk7d5ORMR+++2X3o9feumlbT4Q9Pnnn09//fWvf32Htr2jKm6/Yl1V2bhxY7qvmjZtGgMGDKjX2uriyCOPjG984xsREfHqq6+mZ6+ri379+sWLL76YfgDom2++Gccff3ytLjX96KOPYvDgwem/FT/60Y/il7/8ZZ1rALJHuAFqZVtBoKIFCxak7+2I2DxSUFHFsFPxcpiqVLw3Zs6cOdWu97e//S3efffdWtf45Zdf1jh96wcffBDPPvtsRGweicpVuCnvq23105ZOOeWU9OU4f/rTn9KXpDVt2jTOPffczBZZC+XT8aZSqbjmmmu2u52hQ4emL3m87777ajxAvfPOO7d7OxGb++rEE0+MiIji4uL0dNZV+eCDD+If//hHRER06dIlPS15rhx++OHp0YoJEybErFmzql33vvvui5KSkojY3L8VT0zkg2uvvTb99Y033rhdbfTv3z9efPHF6NixY0REvP7663HCCSfUuP/MmzcvBg8enJ6c4bLLLotbb711u7YPZI9wA9TKQQcdFBdeeGGlqXir8vHHH8fpp5+evqH/0EMP3epSsoqzXE2dOrXG9iqegb711lurnPlq8uTJ8b3vfW+bP8OWfvWrX8WECRO2er24uDjOPvvs9Jn6K664Iho3blzn9jOhvK+WL18eCxcurPX7mjRpku6TKVOmpD+3008/Pdq1a5f5Qrfhkksuid122y0iIv7yl7/EZZddFuvXr692/VWrVsVvf/vbrUYdunbtGmeddVZEbJ7V6/zzz6/yuSy/+MUv4oUXXtjhuq+66qr0/Us//OEPq5zieenSpXHmmWem95cf/vCHOQ8IRUVF6ZvdN27cGGeeeWb6uTcVvf7663HVVVdFxOZRwquvvjqrddbGMcccE//2b/8WEZsvK3vmmWe2q5299947XnjhhfRsda+++moMHTq0yhMHH3/8cQwePDgWLVoUEREXX3xx/Pa3v93OnwDIJlNBA7WyYcOGuOeee+Kee+6JXr16xVFHHRX77bdfdOzYMRo1ahTLli2LN954I5544on07GgtW7aMu+++e6u22rZtG/vvv3+88847MXHixBg5cmQcffTRlW5kLn9g36GHHhoHHXRQTJkyJebPnx99+/aNkSNHRp8+feKLL76IiRMnxl/+8pdIpVIxYsSI+POf/1yrn2fQoEExbdq0OOmkk+LMM8+M4447Lpo3bx7Tp0+PcePGpa+v//rXv57TGZGOPvroePLJJyMi4rTTTosf/OAH0bVr1/Tlal/72teqfNZQxObLz8aMGVPpWUPlz+fItpYtW8YTTzwRgwYNitWrV8edd94Z48ePj7POOisGDBgQu+66a6xevTrmzZsXkydPjhdffDE2bNhQ5fTft912Wzz33HNRXFwcjz76aAwYMCDOP//82GOPPaK4uDj++te/xqRJk+LQQw+NhQsXxuLFi7d7hrLDDz88fvSjH8Wtt94apaWlccQRR8S3v/3tGDRoUBQVFcW0adNi3Lhx6Xt7jjjiiO26N6Q+XH311fH3v/893njjjZgxY0b0798/LrzwwhgwYECsX78+XnrppXjooYfSoezHP/5xHHLIITmuumrXXnttesrlG2+8sdYP9NzS1772tXjhhRdi8ODBsXz58nj55Zdj6NCh8dRTT6VHidesWRODBw+O+fPnR8TmUDRkyJB44oknamy7Q4cOlWaKBHIkBVALxxxzTKqgoCAVEbX6t/fee6emTJlSbXtPPfVUqnHjxtW+v6I5c+akunfvXu26TZs2TY0bNy513333pV+77777ttrmvHnz0su/853vpJ599tlU69atq233oIMOSpWUlFT7M5TX1L179yqX/+xnP0u3NXHixG3WU5XVq1enevfuXW2NVf2cFZ1wwgnpdffaa68a162rbf38Vfnggw9SAwYMqNU+1LRp09TTTz9dZTvvvPNOqmvXrjXufwsWLEh99atfTUVEat99962ynW19RqlUKrVp06bUtddem2rUqFGN9Q4dOjS1evXqan/28vUGDhy4zX6qy7o1WbVqVer444+vse7GjRunrrvuuhrb2Z7PuibV/a5X58svv0z16dMn/Z4t94uKv0u1qXHatGmpdu3apd8zePDg1Oeff75VW3X5t6OfFZAZLksDauWf//xnLFy4MO6555747ne/G1//+tejY8eOUVRUFE2aNIl27drF/vvvH9/73vfiySefjGnTplW6+XtLJ5xwQrz22msxYsSI2GOPPWp8Vsqee+4ZU6dOjZ/85CfRv3//aNasWbRs2TL69OkTl1xySUydOjUuuOCCOv9Mxx57bEybNi2uuOKK6N27dzRv3jxat24dhx56aNx1113x+uuvR/v27evcbia1bNky3nzzzRg1alQccMAB0bp16zqNQgwZMiT9dbYnEqhKv3794p133oknnngizjvvvNhrr71i1113jcaNG0ebNm1iwIABcd5558X9998fn3zySbVn6Pfbb7/44IMP4rrrrou99947/dkdcMABccstt8Rbb70Vu+22W/oyxh25FK+goCBuvvnmmDZtWlx88cXRt2/faNmyZTRr1iy6d+8eZ599dvzjH/+ICRMm5N2U4bvuums8/fTTMWHChBg+fHjsvvvu0bRp02jZsmX07ds3Lr744nj33Xdj9OjRuS61Ro0aNYof//jH6e+3996bcgMGDIjnn38+2rZtGxERL774Ypx00kk7/EwuIPcKUqk6TsEDQINxxBFHxGuvvRZNmjSJjz/+OD1j1M7g/fffj3322SciIi699NIaJ5AAIBmM3AAk1IwZM+K1116LiIhTTz11pwo2EVHpfq9BgwblrhAAska4AUion/70p+mvL7vsshxWknlvvPFGbNy4sdrl99xzT4wdOzYiNs+w9s1vfjNbpQGQQ2ZLA0iIjz76KD766KMoLS2NRx99ND2705AhQ+Lwww/PbXEZdv3118e7774bJ554Yhx44IHRpUuXKCsri7lz58aTTz4ZkydPTq/7+9//Ppo0aZLDagHIlry652bChAkxc+bMWLhwYaxatSo2btwYbdq0if79+8fJJ5+cfkYCAFu74YYbtrrRul27djF58uTo2bNnjqqqH0OGDNnmc2yaNWsWv//97+M73/lOlqoCINfyauTm8ccfj3Xr1kX37t1j9913j4iIRYsWxcsvvxyvv/56XH311bH//vvnuEqA/NaoUaP46le/GkcddVT8/Oc/3+ohqklw1113xYQJE+L555+PuXPnRnFxcaxevTratGkTvXr1iiFDhsTFF18cXbt2zXWpAGRRXo3czJw5M/bcc8+tnuz8z3/+M8aNGxdt27aNsWPHbvfD2AAAgOTKq5TQt2/frYJNxOZnUXTp0iVWrFgRS5YsyUFlAABAvsurcFOT8tGawsK8upIOAADIEw0i3EyaNCmWLFkSXbt23eme0wAAANROXg6DPPnkk7Fo0aJYv359LF68OBYtWhRt27aNyy+/3P02AABAlfIy3Lz77rsxY8aM9Pft27ePSy+9tNYz/lx55ZVbvdaiRYsYPXp0xmoEAADyS16Gm+uvvz4iItauXRsLFy6M8ePHxw033BBnn312nHbaaTmuDgAAyEd5NRV0dcrKyuK6666LefPmxS9+8Yvo1avXDrW3dOnSaN++fURElJSUZKJEtqFDhw4Rob+zSZ9nl/7OLv2dffo8u/R39unz7Crv7yZNmmS03bwcudlSYWFhHHbYYTF37tx4++23dzjcVMxzDSDbJYr+zj59nl36O7v0d/bp8+zS39mnzxu2BnN3fqtWrSIiorS0NMeVAAAA+ajBhJsPPvggIiI6d+6c40oAAIB8lDfh5sMPP4zXX389vvzyy0qvl5WVxdNPPx0vv/xyFBUVxWGHHZajCgEAgHyWN/fcLFu2LO6+++5o1apV7LnnntGqVatYvXp1LFy4MFasWBFNmjSJiy++OH3zEQAAQEV5E2769+8fp556anzwwQexcOHCKC0tjcLCwujUqVMcfPDBceKJJ0aXLl1yXSYAAJCn8ibcdOrUKc4555xclwEAADRQeXPPDQAAwI4QbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQQbgAAgEQozHUB5davXx/vvvtuvP322zFnzpwoLi6OTZs2RZcuXeLggw+OYcOGRbNmzXJdJgAAkKfyZuTm1VdfjVtvvTUmTpwYqVQqBgwYEH379o1PP/00Hnnkkbj22mtj1apVuS4TAADIU3kzclNYWBjHHntsDB06NLp27Zp+fcWKFTFmzJiYN29e3H///XH55ZfnsEoAACBf5c3IzcCBA+PCCy+sFGwiItq2bRsXXHBBRERMnjw5ysrKclEeAACQ5/Im3NSke/fuERGxcePGWL16dY6rAQAA8lGDCDfLli2LiIjGjRtHy5Ytc1wNAACQj/LmnpuaPPXUUxERsd9++0WTJk22uf6VV1651WtFRUUxZsyYiIjo0KFDFBZu/tE7duyYwUqpjv7OPn2eXfo7u/R39unz7NLf2afPs6u8vzMt70dupk6dGhMnTozGjRvH8OHDc10OAACQpwpSqVQq10VU5+OPP47rr78+1q5dG+eff36ceOKJGWn3k08+iQ4dOkRERHFxcUbapGblZ0H0d/bo8+zS39mlv7NPn2eX/s4+fZ5d5f1dm6uy6iJvR26WL18eN998c6xduzaGDRuWsWADAAAkU16Gm9LS0rjpppuipKQkBg0aFOeee26uSwIAAPJc3oWbL774In75y1/G4sWL4+tf/3qMHDkyCgoKcl0WAACQ5/Iq3GzcuDH+8z//M+bMmRMDBgyIK664Iho1yqsSAQCAPJU3yWHTpk3x29/+Nt5///3o169fXHXVVfU2RRwAAJA8eZMennnmmZg8eXJERLRq1SrGjRtX5Xrnnntu7LrrrtksDQAAaADyJtysWbMm/XV5yKnKmWeeKdwAAABbyZtwc9ZZZ8VZZ52V6zIAAIAGKm/uuQEAANgRwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIhbkuACBfHTH+4Xpt/9UzzqnX9gFgZyPcAET9B5nablPgAYDtJ9wAO6VchJna2LIuYQcAak+4AXYK+RpmtkXYAYDaE26AxGqogaYmFX8mQQcAKhNugMRJYqipSvnPKeQAwGbCDZAIO0ugqYrRHADYTLgBGrSdOdRUxWgOADsz4QZokISamgk5AOyMhBugQRFq6kbIAWBnItwADYJQs2OEHAB2Bo1yXQBATY4Y/7Bgk0H6EoAkE26AvOVAvH4IjAAklcvSgLzjwDs7XKoGQNIYuQHyimCTffocgKQQboC84FKp3NL3ACSBcAPknAPr/CBgAtDQCTdATjmYzj8+EwAaKuEGyBkH0fnLZwNAQ2S2NCDrHDg3DGZTA6ChMXIDZJVg0/D4zABoKIQbIGscJDdcPjsAGgLhBsgKB8cNn88QgHwn3AD1zkFxcvgsAchnwg1QrxwMJ4/PFIB8JdwA9cZBcHL5bAHIR8INUC8c/CafzxiAfCPcABnnoHfn4bMGIJ8IN0BGOdjd+fjMAcgXwg2QMQ5yd14+ewDygXADZISDW+wDAOSacAPsMAe1AEA+EG4AyBhBF4BcEm6AHeJgli3ZJwDIFeEG2G4OYqmOfQOAXBBugO3i4BUAyDfCDQD1QgAGINuEG6DOHLRSW/YVALJJuAHqxMEqdWWfASBbhBsAACARhBug1pyBZ3vZdwDIBuEGqBUHp+wo+xAA9U24AQAAEkG4AbbJGXcyxb4EQH0SboAaORgFABoK4QaArBKYAagvwg1QLQeh1Bf7FgD1QbgBAAASQbgBquTMOvXNPgZApgk3AABAIgg3wFacUSdb7GsAZJJwAwAAJIJwA1TiTDrZZp8DIFOEGwAAIBGEGyDNGXRyxb4HQCYINwAAQCIIN0BEOHMOADR8wg0AeUHABmBHCTcAAEAiCDeAM+bkDfsiADuiMNcFVDR37tyYPn16fPTRRzF79uxYsWJFNGnSJB566KFclwYAAOS5vAo348ePj3/961+5LgN2Ks6UAwBJkVfhpnfv3tGjR4/o2bNn9OzZMy666KJclwRAlh0x/uF49Yxzcl0GAA1QXoWbU045JdclAAAADZQJBWAn5pI08pV9E4DtIdwAAACJINwAAACJUJBKpVK5LqI6Z5111nZNBX3llVdu9VpRUVGMGTMmIiI2btwYhYWbbzcqKyvb8ULZJv2dfdvq8z5j78hmObBdZv3gsipf9zcl+/R5dunv7NPn2VXe3wUFBRlt18gNAACQCHk1W1qm3H777TUuLykpiQ4dOkRERHFxcTZK2ul17NgxIvR3NulzkqC6/df+nX36PLv0d/bp8+wq7+8mTZpktF0jN7ATMhMVAJBEwg0AeUsQB6AuhBsAACARhBsAACAR8mpCgalTp8ajjz5a6bWysrIYNWpU+vvTTz89DjjggGyXBonhMh8amiPGPxyvnnFOrssAoAHIq3BTWloas2fPrvRaKpWq9FppaWm2ywIAABqAvAo3gwYNikGDBuW6DAAAoAFyzw0AAJAIwg3sRNxvQ0Nl3wWgNoQbAAAgEYQbAAAgEYQbAAAgEYQb2Em4ZwEASDrhBoAGQUAHYFuEGwAAIBGEGwAAIBGEGwAAIBGEG9gJuFcBANgZCDcANBiCOgA1EW4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4g4cwuBQDsLIQbABoUgR2A6gg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3AABAIgg3kGB9xt6R6xIAALJGuAGgwRHcAaiKcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMAACSCcAMJ1WfsHbkuAQAgq4QbAAAgEYQbAAAgEYQbABokl14CsCXhBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASAThBgAASITCXBewpQ0bNsQTTzwRr732WpSUlETLli1jwIABMXz48Gjfvn2uywMAAPJUXo3cbNiwIUaPHh3jx4+PdevWxYEHHhjt27ePl156KX784x/H0qVLc10iAACQp/Jq5Obxxx+PWbNmRe/eveO6666LZs2aRUTEhAkT4oEHHoixY8fGjTfemOMqAQCAfJQ3IzdlZWXxzDPPRETEBRdckA42ERHDhg2L7t27x4cffhhz587NVYkAAEAey5twM3PmzFi7dm107tw59thjj62WH3zwwRER8a9//SvbpQEAAA1A3oSbBQsWRERUGWwiIvbcc89K6wEAAFSUN/fclJSURERUOyNau3btKq1XkyuvvHKr14qKimLMmDEREdGhQ4coLNz8o3fs2HG76qVu9DdQH/xNyR5/x7NLf2efPs+u8v7OtLwZuVm3bl1ERDRt2rTK5eX34JSvBwAAUFHejNykUqkdWl7R7bffXuPykpKS6NChQ0REFBcX17pdtl/5WRD9DWSSvynZ4+94dunv7NPn2VXe302aNMlou3kzcrPLLrtERMT69eurXF7+esVZ1AAAAMrlTbgpH0lZvnx5lcs/++yzSusBAABUlDfhpnv37hERMW/evCqXlz/fpnw9AACAivIm3PTt2zeaN28ey5YtqzLgvPXWWxERccABB2S7NAAAoAHIm3BTWFgYxx9/fERE3HvvvZVmRZswYUIsWLAg+vbtG7169cpViQAAQB7Lm9nSIiJOO+20mDFjRsyaNSsuv/zy6Nu3b5SUlMTs2bOjVatWcfHFF+e6RAAAIE/lzchNxOYHbf7sZz+L008/PYqKimLKlCnx6aefxsCBA+OWW26JLl265LpEAAAgT+XVyE3E5oAzfPjwGD58eK5LAQAAGpC8GrkBAADYXsINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAACQCMINAA3SrB9clusSAMgzwg0AAJAIwg0AAJAIwg0AAJAIwg0klPsRAICdjXADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADAAAkgnADQIMz6weX5boEAPKQcAMAACSCcAMJ5uw2ALAzEW4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AaFBePeOcXJcAQJ4SbgAAgEQQbiDhnOUGAHYWwg0AAJAIwg0AAJAIwg0AAJAIwg0AAJAIwg0ADYYJMgCoiXADOwEHhADAzkC4AQAAEkG4AQAAEkG4AQAAEkG4AQAAEkG4gZ2ESQVo6OzDAGyLcAMAACSCcAMAACSCcAMAACSCcANA3nO/DQC1IdzATsQBIgCQZMINAACQCIW5LiAiYt26dTF58uT46KOPYvbs2bFgwYIoKyuLESNGxCmnnJLr8gAAgAYgL8LN0qVL43e/+12uywAgD7mcEoDayotw06xZsxg8eHD06tUrevbsGW+99VY89thjuS4LEunVM86JI8Y/nOsyAAAyLi/CTZcuXWLkyJHp76dMmZLDagAAgIbIhAIAAEAiCDewE3IPAw2FfRWAuhBuAACARMiLe24y7corr9zqtaKiohgzZkxERHTo0CEKCzf/6B07dsxqbTsr/Z19+pwkqG7/tX9nnz7PLv2dffo8u8r7O+PtZqKR2267LRYtWlSn91xyySXRq1evTGwe2A6zfnBZ9Bl7R67LgGrN+sFluS4BgAYmI+GmuLg4lixZUqf3rF+/PhObrtLtt99e4/KSkpLo0KFDRGyunfpXfhZEf2ePPqehq2nftX9nnz7PLv2dffo8u8r7u0mTJhltNyPhpvxyLwAAgFwxoQDsxMxERb6ybwKwPYQbAAAgEYQbAAAgEfJmKuhf/epXsXLlyoiIWL58eUREPPvsszFlypSIiGjTpk1cffXVuSoPEuvVM86JI8Y/nOsyIM0laQBsr7wJN/Pnz99qdorly5eng445xwEAgJrkTbi56667cl0CADlm1AaAHeGeG8ABJQCQCMINAACQCMINEBFGb8g9+yAAO0q4AQAAEkG4AdKcOSdX7HsAZIJwAwAAJIJwA1TiDDrZZp8DIFOEGwAAIBGEG2ArzqSTLfY1ADJJuAEAABJBuAGq5Iw69c0+BkCmCTcAAEAiCDdAtZxZp77YtwCoD8INAACQCMINUCNn2Mk0+xQA9UW4ASBrBBsA6pNwA2yTA1IAoCEQboBaEXDYUfYhAOqbcAMAACSCcAPUmjPvbC/7DgDZINwAdeIglbqyzwCQLcINAACQCMINUGfOxFNb9hUAskm4AbaLg1a2xT4CQLYJN8B2c/BKdewbAOSCcAMAACSCcAPsEGfo2ZJ9AoBcEW6AHeZglnL2BQBySbgBMsJBLfYBAHJNuAEyxsHtzstnD0A+EG6AjHKQu/PxmQOQL4QbALabYANAPhFugIxzwLtz8DkDkG+EG6BeOPBNNp8vAPlIuAHqjQPgZPK5ApCvhBugXjkQThafJwD5TLgB6t2rZ5zjoDgBfIYA5DvhBsgaB8cNl88OgIZAuAGyykFyw+MzA6ChEG6ArHOw3HD4rABoSIQbICccNOc/nxEADU1hrgsAdl7lB89HjH84x5VQkVADQENl5AbIOQfT+cNnAUBDJtwAecFBde75DABo6FyWBuQNl6nlhlADQFIYuQHyjoPt7NHXACSJkRsgLxnFqV9CDQBJZOQGyGsOwjPr1TPO0acAJJZwA+Q9B+SZoQ8BSDqXpQENhkvVto9QA8DOQrgBGhwhp3aEGgB2NsIN0GAJOVUTagDYWQk3QINX8WB+Zw06Ag0ACDdAwuxsozlCDQD8L+EGSKQkj+YINABQNeEGSLwkBB2BBgC2TbgBdipbhoR8DTvCDADUnXAD7NSqChHZDjyCDABkhnADsIVthY26hh/hBQCyQ7gBqKOqwkrHjh0jIqK4uDjb5QAA/6NRrgsAAADIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIBOEGAABIhMJcFxARsXjx4pgyZUpMnz49Pvnkk1i1alW0aNEi+vTpE0OHDo1+/frlukQAACDP5UW4GT16dHz22Wexyy67xF577RW9e/eOjz/+OCZPnhxTpkyJ8847L4YOHZrrMgEAgDyWF+GmW7du8e1vfzsOOeSQKCz835Kee+65+NOf/hQPPvhgDBgwILp165bDKgEAgHyWF/fcXHfddXHEEUdUCjYREcccc0wMGDAgNm3aFG+88UaOqgMAABqCvAg3NenevXtERKxYsSLHlQAAAPks78PNsmXLIiKiTZs2uS0EAADIa3kdbpYuXRpTp06NiIgDDzwwx9UAAAD5rCCVSqVyXURVvvzyy7jxxhtj5syZcdhhh8UVV1xR6/deeeWVW71WVFQUY8aMiYiIjRs3pu/vKSsry0i91Ex/Z58+zy79nV36O/v0eXbp7+zT59lV3t8FBQWZbTcTjdx2222xaNGiOr3nkksuiV69elW7/N57742ZM2dG586d48ILL9zREgEAgITLSLgpLi6OJUuW1Ok969evr3bZ+PHj47nnnovWrVvHqFGjomXLlnVq+/bbb69xeUlJSXTo0CEiNtdO/evYsWNE6O9s0ufZpb+zS39nnz7PLv2dffo8u8r7u0mTJhltNyPhpvxyr0x45pln4pFHHonmzZvHqFGjokuXLhlrGwAASK68mlDglVdeifvuuy+aNm0a11xzTfTo0SPXJQEAAA1E3oSbqVOnxt133x2NGzeOq666Kvr27ZvrkgAAgAYkL8LNzJkz0/fJXHHFFTFgwIAcVwQAADQ0GbnnZkfdcsstsWHDhujUqVNMmTIlpkyZstU6ffv2jaOPPjoH1QEAAA1BXoSbtWvXRkTEp59+Gp9++mm16wk3AABAdfIi3DzyyCO5LgEAAGjg8uKeGwAAgB0l3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIkg3AAAAIlQmOsCIiIWLFgQzz//fMydOzdKSkpi9erV0aRJk+jWrVsceeSRccwxx0Tjxo1zXSYAAJDH8iLcfPjhh/Hss89Gx44do1u3btGqVatYvXp1zJw5M2bPnh2TJ0+On/zkJ1FYmBflAgAAeSgv0sL+++8fd955Z3Tu3LnS6ytXrozRo0fHe++9Fy+++GIce+yxOaoQAADId3lxz03nzp23CjYREW3atIlTTjklIiLee++9LFcFAAA0JHkRbmrSqNHmEl2SBgAA1CSvE8OaNWtiwoQJEbH50rVMKSgoqPJr6p/+zj59nl36O7v0d/bp8+zS39mnzxu2glQqlcp1EeU++eSTeOyxxyKVSsWqVati1qxZsW7duhgyZEh8//vfr/XOduWVV271WosWLWL06NGZLhkAAMgTeTVys2rVqpg0aVKl144//vg4++yzpWgAAKBGGRm5ue2222LRokV1es8ll1wSvXr1qnLZpk2boqSkJCZPnhx//etfo02bNjFq1Kjo1KnTjpaads0110RExJgxYzLWJtXT39mnz7NLf2eX/s4+fZ5d+jv79Hl21Vd/Z2Tkpri4OJYsWVKn96xfv77aZY0aNYpOnTrFsGHDolOnTnHrrbfGvffem+6ETNiwYUPG2mLb9Hf26fPs0t/Zpb+zT59nl/7OPn2eXfXV3xkJN/WZcA866KBo1qxZTJs2LcrKysyaBgAAVCnvp4IuKCiIli1bxqZNm2LNmjW5LgcAAMhTeR9uli1bFsuXL49ddtkldt1111yXAwAA5Km8CDd/+9vfYtmyZVu9vmTJkrjjjjsilUrFwIED0w/0BAAA2FJePOfmP/7jP6KkpCR69OgRnTt3jojNkxTMnTs3UqlU9OvXL6699tpo1qxZjisFAADyVV6Em1deeSXeeeedmDNnTqxcuTI2bNgQLVu2jB49esThhx8eRx11lFEbAACgRnkRbgAAAHaU4RAAACARhBsAACARhBsAACARhBsAACARhBsAACARCnNdQD4bP358PPLIIxERcfnll8fhhx+e44qSY8GCBfH888/H3Llzo6SkJFavXh1NmjSJbt26xZFHHhnHHHNMNG7cONdlJsbixYtjypQpMX369Pjkk09i1apV0aJFi+jTp08MHTo0+vXrl+sSE2fdunUxefLk+Oijj2L27NmxYMGCKCsrixEjRsQpp5yS6/IarA0bNsQTTzwRr732WpSUlETLli1jwIABMXz48Gjfvn2uy0uUuXPnxvTp09P78IoVK6JJkybx0EMP5bq0RFq/fn28++678fbbb8ecOXOiuLg4Nm3aFF26dImDDz44hg0b5nl/9WDChAkxc+bMWLhwYaxatSo2btwYbdq0if79+8fJJ58cu+22W65LTKw1a9bEFVdcEaWlpfGVr3wlfvOb32SkXeGmGkuWLInHH388CgoKwmzZmffhhx/Gs88+Gx07doxu3bpFq1atYvXq1TFz5syYPXt2TJ48OX7yk59EYaFdNBNGjx4dn332Weyyyy6x1157Re/evePjjz+OyZMnx5QpU+K8886LoUOH5rrMRFm6dGn87ne/y3UZibJhw4YYPXp0zJo1K9q2bRsHHnhgFBcXx0svvRRTp06Nm266Kbp06ZLrMhNj/Pjx8a9//SvXZew0Xn311fjDH/4QERG77bZbDBgwIL744ov4f//v/8UjjzwSr732Wtxwww3RunXrHFeaLI8//nisW7cuunfvHrvvvntERCxatChefvnleP311+Pqq6+O/fffP8dVJtN//dd/xerVqzPeriPHKqRSqfjDH/4QzZs3j7322ssf93qw//77x5133hmdO3eu9PrKlStj9OjR8d5778WLL74Yxx57bI4qTJZu3brFt7/97TjkkEMqBcbnnnsu/vSnP8WDDz4YAwYMiG7duuWwymRp1qxZDB48OHr16hU9e/aMt956Kx577LFcl9WgPf744zFr1qzo3bt3XHfddemz2BMmTIgHHnggxo4dGzfeeGOOq0yO3r17R48ePaJnz57Rs2fPuOiii3JdUqIVFhbGscceG0OHDo2uXbumX1+xYkWMGTMm5s2bF/fff39cfvnlOawyea6++urYc889o6ioqNLr//znP2PcuHHx+9//PsaOHeth8hk2Y8aMmDRpUgwZMiSef/75jLbtk6rCCy+8EB9++GGcd9550aJFi1yXk0idO3feKthERLRp0yZ9yc57772X5aqS67rrrosjjjhiq5GwY445JgYMGBCbNm2KN954I0fVJVOXLl1i5MiRMWTIkNhjjz38x7iDysrK4plnnomIiAsuuKDS5TnDhg2L7t27x4cffhhz587NVYmJc8opp8RZZ50V//Zv/xZt2rTJdTmJN3DgwLjwwgsrBZuIiLZt28YFF1wQERGTJ0+OsrKyXJSXWH379t0q2EREHHvssdGlS5dYsWJFLFmyJAeVJdeGDRviT3/6U3Tr1i2++c1vZrx9/9tuYeXKlfHQQw/F1772tTjyyCNzXc5Oqfwg0CVp2dG9e/eI2Hx2EPLVzJkzY+3atdG5c+fYY489tlp+8MEHR0QYaSeRyv9Ob9y4sV4u46Fqjkfqx1//+tdYtmxZXHjhhfVyf7Vws4V77703NmzYEBdeeGGuS9kprVmzJiZMmBAR4RrXLFm2bFlEhDOz5LUFCxZERFQZbCIi9txzz0rrQZKU/51u3LhxtGzZMsfV7BwmTZoUS5Ysia5du0anTp1yXU5iLFiwICZMmBCDBg2K/v3718s2RNEK3n777XjzzTfjrLPO2mpYmPrxySefxGOPPRapVCpWrVoVs2bNinXr1sWQIUPiiCOOyHV5ibd06dKYOnVqREQceOCBOa4GqldSUhIRUe2MaO3atau0HiTJU089FRER++23XzRp0iTH1STTk08+GYsWLYr169fH4sWLY9GiRdG2bdu4/PLLXVacIZs2bUrf0/7tb3+73rYj3PyPdevWxbhx46Jr165x8skn57qcncaqVati0qRJlV47/vjj4+yzz46CgoIcVbVz+PLLL+Puu++OjRs3xmGHHZY+8w35aN26dRER0bRp0yqXl9+DU74eJMXUqVNj4sSJ0bhx4xg+fHiuy0msd999N2bMmJH+vn379nHppZf6vzGDnnnmmfjoo4/i4osvjlatWtXbdhITbm677bZYtGhRnd5zySWXRK9evSIi4s9//nMsX748fvrTnzorUgs72t/l+vbtG4888khs2rQpSkpKYvLkyfHXv/413n333Rg1apSh4P+Rqf6u6N57742ZM2dG586dXYZZhfroc7bftqbkN2U/SfTxxx/HnXfeGalUKs4999zo0aNHrktKrOuvvz4iItauXRsLFy6M8ePHxw033BBnn312nHbaaTmuruErKSmJv/zlL9G/f/8YNGhQvW4rMeGmuLi4zrNZrF+/PiIiPvroo3j22WfjqKOOin322ac+ykucHenvqjRq1Cg6deoUw4YNi06dOsWtt94a9957b1xzzTU7WmoiZLq/x48fH88991y0bt06Ro0a5RruKmS6z9kxu+yyS0RU38flr3vIIUmxfPnyuPnmm2Pt2rUxbNiwOPHEE3Nd0k6hRYsW0a9fv7j22mvjuuuui//+7/+Offfd14mrHTRu3LgoKyvLysnUxISbMWPGbPd7p06dGqlUKhYuXBg33HBDpWWLFy+OiP89GDzkkEPi+OOP35FSE2FH+ntbDjrooGjWrFlMmzYtysrKzFISme3vZ555Jh555JFo3rx5jBo1ykMPq1Gf+zh116FDh4jYfMBXlc8++6zSetCQlZaWxk033RQlJSUxaNCgOPfcc3Nd0k6nsLAwDjvssJg7d268/fbbws0Omjp1arRo0SLGjRtX6fWNGzdGxOaRnfJj8GuuuWaHTlQ5aqxg/vz51S5bvHhxLF682JBwFhQUFETLli2jpKQk1qxZYxavDHrllVfivvvui6ZNm8Y111xjf6bBKJ8Kd968eVUuL3++Tfl60FB98cUX8ctf/jIWL14cX//612PkyJHuQc2R8vtCSktLc1xJMqxduzY++OCDKpdt2LAhvezLL7/coe0INxFx1llnxVlnnVXlsrvuuismTZoUl19+eRx++OFZrmzntGzZsli+fHnssssuseuuu+a6nMSYOnVq3H333dG4ceO46qqrom/fvrkuCWqtb9++0bx581i2bFnMmzdvqymh33rrrYiIOOCAA3JRHmTExo0b4z//8z9jzpw5MWDAgLjiiivM1JVD5QfbVT10nLp55JFHqnz9008/jUsuuSS+8pWvxG9+85uMbMtvDDnxt7/9LT1vf0VLliyJO+64I1KpVAwcONAf9QyZOXNm3H777RERccUVV8SAAQNyXBHUTWFhYfqS4HvvvbfSrGgTJkyIBQsWRN++fV06QoO1adOm+O1vfxvvv/9+9OvXL6666iqXZdezDz/8MF5//fWtRgrKysri6aefjpdffjmKiorisMMOy1GFbA+/NeTEP//5z/jzn/8cPXr0SJ8RKS4ujrlz50YqlYp+/frFiBEjclxlctxyyy2xYcOG6NSpU0yZMiWmTJmy1Tp9+/aNo48+OgfVJdevfvWrWLlyZUT8770izz77bLr/27RpE1dffXWuymtwTjvttJgxY0bMmjUrLr/88ujbt2+UlJTE7Nmzo1WrVnHxxRfnusREmTp1ajz66KOVXisrK4tRo0alvz/99NONlmXIM888E5MnT46IzZdDbXlvQrlzzz3XVQ0ZsmzZsrj77rujVatWseeee0arVq1i9erVsXDhwlixYkU0adIkLr74YvfyNTDCDTlx9tlnxzvvvBNz5syJd999NzZs2BAtW7aMfffdNw4//PA46qijjNpk0Nq1ayNi8/Dvp59+Wu16wk1mzZ8/P4qLiyu9tnz58nTQ6dixYy7KarCKioriZz/7WTz++OPx6quvxpQpU6JFixYxcODAGD58uAOQDCstLY3Zs2dXei2VSlV6zb0ImbNmzZr01+UhpypnnnmmcJMh/fv3j1NPPTU++OCDWLhwYZSWlkZhYWF06tQpDj744DjxxBNNutMAFaQ8HAAAAEgAp8YBAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBEEG4AAIBE+P8frHYNxRvqCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Coefficient for Heun's\n",
    "c2, a21 = 1, 1\n",
    "b1, b2 = 0.5, 0.5\n",
    "\n",
    "# Coefficient for Midpoint\n",
    "#c2, a21 = 1/2, 1/2\n",
    "#b1, b2 = 0.0, 1.0\n",
    "\n",
    "k1 = ff(t, y)\n",
    "k2 = ff(t + c2*h, y+ h*a21*k1)\n",
    "\n",
    "# Factor expression\n",
    "expr = sy.simplify((y+ h*(b1*k1 + b2*k2)) / y)\n",
    "print(expr)\n",
    "\n",
    "# Sig\n",
    "z = sy.Symbol('z')\n",
    "sigf = sy.lambdify(z, expr.subs(h*lam, z), 'numpy')\n",
    "\n",
    "# Make grid\n",
    "xx = np.linspace(-3, 3, 201)\n",
    "yy = np.linspace(-3, 3, 201)\n",
    "X, Y = np.meshgrid(xx, yy)\n",
    "z = X + Y*1j\n",
    "\n",
    "# Amplication factor of Explicit Euler (z=lambda h)\n",
    "sig = sigf(z)\n",
    "\n",
    "# Same scale of x and y\n",
    "plt.axis('equal')\n",
    "\n",
    "# Stability region\n",
    "plt.contourf(X,Y,abs(sig), levels=[0, 1])\n",
    "\n",
    "# Title\n",
    "plt.title('Stability region of RK2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06d51a6d",
   "metadata": {},
   "source": [
    "### 4th-order Runge-Kutta Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d8ea4cf7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0416666666666667*h**4*lambda**4 + 0.166666666666667*h**3*lambda**3 + 0.5*h**2*lambda**2 + 1.0*h*lambda + 1.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Stability region of RK4')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzcAAAKOCAYAAACbY+fTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAABcSAAAXEgFnn9JSAABJ/0lEQVR4nO3deXxU1f3/8XfIxhIgbCFYMOwEUCMWRAQhsqmAiopsdatSSykIRfmJgoqFWqxKLSq4IFj9KhZBUCmCIIuCSCIhgGAoWxLWkCCQgJAQuL8/aKYJ2ZOZuXfOvJ6PB49HMvfOmc+cmQz3PefccwMsy7IEAAAAAD6uit0FAAAAAIA7EG4AAAAAGIFwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEGwAAAABGINwAAAAAMALhBgAAAIARCDcA/EZycrICAgIUEBCghx56qNLtNW3aVAEBAWratGmR26dMmeJ6vLVr13q8HpRfaa8R/ufs2bOaPn26OnfurPDwcFWpUsXVdydPnrS7PACQJAXZXQAA3/P999/r//7v/7Rx40YlJycrMzNTISEhqlu3rlq0aKGYmBjdcMMN6t27txo0aFDo/snJyXrvvfckSbGxsYqNjfXuE/Ahr776qk6ePKnw8HCNGzfO7nLgp3755Rd1795dmzdvdkt7AQEBxW4LCwtT/fr1dc0112jAgAEaPny4atSoUWJ7ycnJatasmSQpKipKycnJJe5/8uRJ9evXTxs3bpQktWnTRqtWrVLjxo1LrX3v3r26+uqrdfbsWddtlmWVej8A3kG4AVBmp06d0ogRI7Rw4cJC23Jzc/XLL7/o4MGDWrdunWbOnKmAgABlZWUVOjBJTk7W888/7/qdcFO8V199VSkpKYqKiiLcwDZvvvmmK9hcddVVevTRR3XFFVcoMDBQkkoNH+Vx+vRpnT59WsnJyfr88881bdo0ffzxx+rSpYtb2k9PT1ffvn2VmJgoSbr22mu1YsUKRURElOn+v/vd7woEGwDOQrgBUCbnz5/XLbfcok2bNkmSQkJCdOedd6pr166KjIyUZVk6evSoEhMTtWrVKh06dEiWZRn9jWZp3w6XpmnTpkb3jy+YMmWKpkyZYncZjrds2TJJl0ZcVqxYoSuuuMJtbS9evLjA75mZmdqyZYs++OADHT9+XKmpqerXr58SExMVFRVVqcc6dOiQevfuraSkJElSly5dtGzZMoWHh5fp/nPmzNGaNWtUo0YNnTlzplK1APAMwg2AMnnjjTdcwaZ58+Zavny5WrVqVeS+lmXp+++/1+zZs1WlCqf2Ab7uwIEDkqSGDRu6NdhI0sCBAwvd9sADD+ipp55Sjx49lJSUpJMnT2ratGl65513Kvw4+/btU69evVxfSvTq1UufffZZmUedjhw5ogkTJkiS/vznP+vxxx+vcC0APIejDgBl8uGHH7p+nj17drHBRrr07W6XLl30/vvvq3r16t4oD4AHZWdnS5JCQ0O99pgRERF6+eWXXb9/8cUXFW5r586duummm1zB5o477tC///3vck2n++Mf/6iTJ0+qQ4cOGjt2bIVrAeBZhBsAZZI3jUOSunfvXqE21q5dq4CAAN18882u255//nnXikv5/10uNTVVr7/+uu699161adNGYWFhCgkJUUREhGJjY/Xiiy/q1KlT5a4pOTlZ48ePV5s2bVSjRg3VrVtXXbt21ezZs3XhwoUS71vaamlleeziVkvLazslJUWSlJKSUmQ/5S3M0LlzZwUEBCgkJERpaWmlPnZaWpqCg4MVEBCgG264oUL1X/78z5w5o1deeUVdunRRRESEqlSpUuS38ufPn9fcuXN1xx13qEmTJqpatarCw8N1zTXX6PHHHy/zdL+TJ0/q2Wef1TXXXKOwsDDVqVNHHTt21EsvveSaMlTZFe3y+/HHHzV69Gi1a9dOtWvXVrVq1dSsWTP95je/0ZdffllqvXmPk3eO2ZkzZ/Tyyy+rY8eOCg8PV40aNXTVVVdp0qRJbl99bNmyZRo+fLiaNm2qatWqqXbt2mrfvr3GjBmjHTt2FHmf9957z1VzSe/DvPegJ9x0002un9PS0pSZmVnuNhISEtSjRw8dPnxYkjRs2DAtWrSoXEFt0aJFWrx4sapUqaK3337bda4RAAeyAKAMqlatakmyJFkpKSkVamPNmjWuNkr7d/n9AgICSr1PgwYNrG+//bbYx9+/f79r3wcffNBasWKFVbt27WLb69Spk5WRkVFse1FRUZYkKyoqqsjtzz33nKutNWvWlFpPUW2X9m/evHmWZVnW3LlzXbe9+OKLxdacZ/r06a7958yZU+r+Rcn//Pfu3WtFR0cXqu/OO+8scJ/NmzdbzZs3L/E5hYSEWG+++WaJj52YmGg1atSo2Dbat29vpaamVvo1sizLunjxovX0009bVapUKbHuAQMGWFlZWcXWnLdfjx49rD179lht27Yttq2mTZtaycnJJfZBWWRmZlq33XZbiXUHBgZazzzzTKH7zps3r1zvwfIo7m/9cufOnSuw79GjR4vcL//fUv7Xev369QX+xh999FHrwoUL5ar1xIkTVmRkpCXJGjt2bLmfAwDv4pwbAGXSsmVL/fjjj5Kk1157TS+99FK527jqqqu0ePFi/fjjj3rmmWckSUOGDNHQoUNLvN+5c+dkWZbat2+vm2++WW3btlW9evV07tw5HThwQEuWLNHmzZuVnp6uAQMGKDExsdTRlJSUFA0ZMkSZmZm69957dcstt6h69eratm2b5syZo4yMDMXHx6t///7asGGD17+pffvtt/XLL7/o0UcfVXp6uho0aKC333670H7XXXedJGno0KEaP368Tp48qTlz5uj//b//V2zblmXp3XfflSTVrFlTQ4YMqVSt2dnZuvvuu5WUlKQePXrorrvuUqNGjZSWlqasrCzXfps2bVKvXr1coyq9evXSbbfdpiZNmujcuXPauHGj3n//ff3yyy8aOXKkQkNDi7z+z9GjR9WnTx+lp6dLkqKjo/XQQw+padOmysjI0CeffKJ169ZpyJAhys3NrdRzk6Qnn3zS9X4PCgrS8OHDFRsbq9DQUG3ZskXvvvuuTpw4oaVLl6pfv35as2ZNie+XzMxM9e/fX7t27dIdd9yh2267TXXr1tW+ffs0a9YsHThwQMnJyXrooYe0Zs2aCtedm5urW265xbXccd26dfXII4/o2muvVXZ2ttasWaOPPvpIFy5c0NSpU5Wbm6sXXnjBdf+ePXu6TvYv6X2Y9x70hPyjSqGhoWVe0UySVq5cqYEDB+qXX36RJD3++OMFprmV1eOPP66jR4+qcePGmjZtWrnvD8DL7E5XAHzD1KlTC30j/9lnn1knT54sd1v5R3Cee+65UvdPTk62tm3bVuI+8+fPtwIDAy1J1kMPPVTkPvm/3ZVkBQUFWUuWLCm0X1pamtW+fXvXfq+88kqR7Xly5Kasj5HfmDFjSh2FsCzLWr16dYFvsivq8tGlmTNnFrtvVlaWa/8aNWpYy5YtK3K/3bt3W1deeaUlyapevbp17NixQvsMHz7c9Zj33HOPlZ2dXWifadOmFaitoq/R+vXrXaOGNWvWtDZs2FBon8OHDxcYtZo+fXqRj5W/npCQEOuLL74otE96errVtGlT135xcXFFtlUWf/nLX1zttGvXzjpy5Eihfb755hurRo0aliSrSpUq1nfffVdkW+V5H5ZF/r4oyaBBgwqMeBXn8pGbJUuWWKGhoa7bpkyZUqE6V61a5Wrj8s+Ksj4HAN7FXySAMjlz5ozVsWPHQlNSAgICrFatWlnDhg2zZs6cae3YsaPUtsobbsrqwQcftCRZ1apVs3JycgptvzzcTJw4sdi2tm/f7gpLTZo0sXJzcwvt47Rws2PHDld7v/nNb4rdL384iI+PL7Xd4uQPN4MGDSpx3xkzZrj2ff/990vc9+uvv3btO23atALbDh8+bAUFBVmSrIYNG1qZmZnFttOzZ89Kh5s777zTtf2tt94q9rG2b9/uqisyMrLIwJX/vffnP/+52Lbefvtt135Tp04tdr+SZGdnWxEREa4QX9Lf5ezZs12Pd9dddxW5jzfDTWZmprVu3Trr9ttvL7Dfl19+WWx7+f+WQkNDXa9FQECA9fe//71CNZ45c8Y1hbKofiHcAM7EggIAyqR69epau3atxowZo5CQENftlmVp9+7dmj9/vh577DG1b99eHTp00KJFi7xeY95F/s6ePatt27aVuG9gYGCJF8W86qqrdMstt0i6tAxufHy82+r0lHbt2rlOwF60aJFOnDhRaJ+ff/5Zn376qSSpQ4cO6tixo1see/To0SVuf//99yVJjRo10m9+85sS9+3Zs6drueGvvvqqwLZly5a5ppr99re/Vc2aNYttZ8yYMaXWXZLs7GzXQgH169fXb3/722L3veqqq3TbbbdJujRtLm8qWFECAwNL7K9evXq5fi7uZP/SbNiwQceOHZMk9e/fX+3atSt234cfflj16tWTdKl/c3JyKvSYFXX5AgW1atVSjx49CqyONmPGDN16661lai87O9v1Hrnyyiv1wAMPVKiuZ599Vvv27VPNmjX12muvVagNAN5HuAFQZjVq1NDMmTN18OBBzZ49W/fcc49+9atfFdovMTFRgwYN0sMPP6yLFy+67fE3bdqkMWPGqFOnTqpXr55CQkIKHBSNHDnSte/BgwdLbKt9+/Zq2LBhifv07NnT9bMvhBtJ+v3vfy/p0nlKH3zwQaHtH3zwgc6dOyfp0pXW3SEoKEidO3cudntmZqYrbDZq1Eiff/65lixZUuK/sLAwSQVX6ZOkH374wfVz3qpjxSlte2kSExNdB/qxsbEKDg4ucf8+ffq4fs67JlRRWrdurTp16hS7Pf/fVFEBtSzi4uKKrKsoISEh6tGjh6RLwSAxMbFCj+kJHTp00Pbt2/WnP/2pzPcJDw939WFKSop69+5d7n784Ycf9Oqrr0qSXnjhhSI/5wA4EwsKACi3Bg0aaOTIka4wcezYMW3atElffvmlPvzwQ9dyrfPmzVOLFi00adKkSj1eTk6ORowYUeTBenFKWzK2ZcuWpbaRf5+8ZWSdbtCgQRo7dqyOHz+uOXPm6LHHHiuwPW8hgerVq2v48OFuecy6deuqatWqxW5PTU11hdyEhATdddddZW77559/LvD7kSNHXD83a9asxPuGh4erTp06FQ4I+R+rpOs65WndurXr56NHjxa7X/369UtsJ/8SxXlBtLw8Vbsn5C1aIF0adU1OTtaHH36oHTt2aMuWLXr99dc1a9asMl8QuHbt2lq5cqViY2N1+PBhbdmyRX369NGqVasUHh5e6v3Pnz+vRx55RBcuXND111+vUaNGVfSpAbABIzcAKi0iIkK33367Zs2apb1796pbt26ubS+++KLOnj1bqfb/+Mc/uoJNaGio7rrrLr3wwgv65z//qYULF2rx4sVavHhxgWlIpV2jpiwXF81/gb/Tp09XsHrvyr/C2Pbt2wuMIHz//ffavn27JGnw4MGqXbu2Wx6zWrVqJW6vyPWH8ly+2ln+16EsF2Asz0UaL1eZx8q/StzlynqQXhmeqt0TBg4c6Po3bNgwPfXUU9q+fbvr7/mtt95yra5YVq1atdKaNWvUqFEjSdLmzZvVt2/fMr0X//a3v2nbtm0KCgrS22+/7ZXXC4D7MHIDwK3q16+v+fPnq1mzZsrNzVVWVpbi4uJc017KKyUlxTXa0LhxY61bt07Nmzcvct9Dhw6Vud285WFLkrdksSTXNClf8Pvf/14zZsyQZVl65513XFPG3nnnHdc+7pqSVhb5++6hhx7SvHnzKtxW/oPwsryGZdmnOPnrLu/7paRzgbzBl2uXLp2H8/e//10bN27UDz/8oOnTp+vOO+/U9ddfX+Y2WrdurTVr1ig2NlZHjx5VfHy8brnlFn311VeqVatWsfebM2eOpEsjt1988UWBc3+Kk3+J6AkTJpTrAqEA3ItwA8DtGjdurNatW2vnzp2SKjel6+uvv5ZlWZKkiRMnFhtsJLmuol4We/bsKdc+eSe4+4JWrVqpZ8+e+vrrr/Wvf/1Lr776qizL0oIFCyRdOt/oxhtv9Fo9+fsu71pJFZX3Tbwk7d+/v8QpVydPniw0ra2ij7V79+5S98+/T/772uHy2ks778ZJtecJDAzUjBkz1L17d128eFGPP/64vv3223K10aZNG61evVo333yz0tLStGnTJlfAKS7E5X3eJCUllXnEKP9+o0ePJtwANmKsFYBH5F9R7fJRj/zTPPIOJIqTlpbm+rlFixYl7nv5ylol2bFjR4G2i5L/AoqdOnUqc9vulNdXpfXT5fIWFjh9+rTmz5+v+fPnu6YqeXPURrp0jlbbtm0lXZoeVJ4Qern8q7utXbu2xH3XrVtX4ceRpGuvvdb1Pl67dm2pFwRdtWqV6+fyjDB4Qv7Hz19XUc6fP+/qq9DQUMXExHi0tvK46aabdPPNN0uS1q9f71q9rjzatm2r1atXuy4A+v333+vWW2/1mammAMqHcAOgTEoLAvmlpKS4zu2QLo0U5Jc/7OSfDlOU/OfG7N27t9j9PvvsM23durXMNV64cEEzZ84sdvvOnTu1YsUKSZdGouwKN3l9VVo/XW7gwIGKjIyUdGk6Wt6UtNDQUN1///3uLbIM8pbjtSxLEydOrHA7/fv3V1DQpUkH8+bNK/EAtbLL94aGhqpfv36SpPT0dNdy1kXZuXOn/v3vf0uSIiMjXcuS26Vr166u1QCXLl2qXbt2FbvvvHnzlJGRIelS/+b/YsIJnnrqKdfPzz//fIXaaNeunVavXq0GDRpIkr777jvddtttRb5/kpOTZV26DmCJ//LLf3tZFi0A4DmEGwBl0qlTJ40YMaLAUrxFOXjwoO655x7XCf1dunQpNJUs/ypXCQkJJbaX/xvol19+uciVr+Li4vTwww+X+hwu99JLL2np0qWFbk9PT9fQoUNd39SPGzdOgYGB5W7fHfL66vjx40pNTS3z/YKDg119Eh8f73rd7rnnHtWtW9f9hZZi9OjRatKkiSTp448/1mOPPabs7Oxi9z916pT+8Y9/FBp1aNSokQYPHizp0qpeDz30UJHXZfnLX/6ir7/+utJ1P/HEEwoICJAk/elPfypyieejR4/q3nvvdb1f/vSnP9keEEJCQjR27FhJl0Zm7r33Xtd1b/L77rvv9MQTT0i6NEo4YcIEr9ZZFn369NGvf/1rSZeW2F6+fHmF2mnfvr2+/vpr12p169evV//+/cv9xQEAZ+OcGwBlkpOTo3fffVfvvvuuWrZsqe7du+vaa69VgwYNVKVKFaWlpWnjxo1asmSJa3W0sLAwzZo1q1BbderUUYcOHbRlyxatWbNGI0eOVK9evQrMgc+7YF+XLl3UqVMnxcfHKzk5WdHR0Ro5cqTatGmjs2fPas2aNfr4449lWZaGDx+ujz76qEzPJzY2VomJibrjjjt077336pZbblH16tW1bds2zZkzR+np6ZIuhau8g0Q79OrVS59//rkk6e6779Yf/vAHNWrUyDVd7eqrry72Ghy/+93vNH369ALXGnr00Uc9X3QRwsLCtGTJEsXGxiorK0uvvfaaFi5cqMGDBysmJka1atVSVlaW9u/fr7i4OK1evVo5OTlFLv/9yiuvaOXKlUpPT9eiRYsUExOjhx56SM2aNVN6ero++eQTrVu3Tl26dFFqaqoOHTpU4RWvunbtqscff1wvv/yyMjMz1a1bN913332KjY1VSEiIEhMTNWfOHNe5Pd26ddPjjz9eqb5ylwkTJuiLL77Qxo0btX37drVr104jRoxQTEyMsrOztXbtWn344YeuUPbkk0/qhhtusLnqoj311FMaNGiQpEujN2W9oOflrr76an399dfq2bOnjh8/rm+++Ub9+/fXsmXLyrSCIgAfYAFAGfTp08cKCAiwJJXpX/v27a34+Phi21u2bJkVGBhY7P3z27t3rxUVFVXsvqGhodacOXOsefPmuW6bN29eocfcv3+/a/uDDz5orVixwqpdu3ax7Xbq1MnKyMgo9jnk1RQVFVXk9ueee87V1po1a0qtpyhZWVlW69ati62xqOeZ32233ebat1WrViXuW16lPf+i7Ny504qJiSnTeyg0NNT68ssvi2xny5YtVqNGjUp8/6WkpFi/+tWvLEnWNddcU2Q7pb1GlmVZFy9etJ566imrSpUqJdbbv39/Kysrq9jnnrdfjx49Su2n8uxbklOnTlm33npriXUHBgZakydPLrGdirzWJSnub704Fy5csNq0aeO6z+Xvi/x/S2WpMTEx0apbt67rPj179rR++eUXjz4HAN7BtDQAZfLVV18pNTVV7777rn7729/q+uuvV4MGDRQSEqLg4GDVrVtXHTp00MMPP6zPP/9ciYmJBU7+vtxtt92mDRs2aPjw4WrWrFmJ10pp3ry5EhIS9PTTT6tdu3aqWrWqwsLC1KZNG40ePVoJCQl65JFHyv2c+vbtq8TERI0bN06tW7dW9erVVbt2bXXp0kVvvPGGvvvuO9WrV6/c7bpTWFiYvv/+e02aNEnXXXedateuXa5RiN69e7t+9vZCAkVp27attmzZoiVLluiBBx5Qq1atVKtWLQUGBio8PFwxMTF64IEH9N577+nIkSPFfkN/7bXXaufOnZo8ebLat2/veu2uu+46vfjii9q0aZOaNGnimsZYmal4AQEBeuGFF5SYmKhRo0YpOjpaYWFhqlq1qqKiojR06FD9+9//1tKlSx23ZHitWrX05ZdfaunSpRoyZIiuvPJKhYaGKiwsTNHR0Ro1apS2bt2qqVOn2l1qiapUqaInn3zS9XtFz73JExMTo1WrVqlOnTqSpNWrV+uOO+6o9DW5ANgvwLLKuQQPAMBndOvWTRs2bFBwcLAOHjzoWjHKH+zYsUNXXXWVJGnMmDElLiABADADIzcAYKjt27drw4YNkqS77rrLr4KNpALne8XGxtpXCADAawg3AGCoZ5991vXzY489ZmMl7rdx40adP3++2O3vvvuuZs+eLenSCmu33367t0oDANiI1dIAwBB79uzRnj17lJmZqUWLFmnJkiWSLp1307VrV3uLc7NnnnlGW7duVb9+/dSxY0dFRkYqNzdX+/bt0+eff664uDjXvm+++aaCg4NtrBYA4C2OOudm6dKlSkpKUmpqqk6dOqXz588rPDxc7dq105133um6RgIAoLApU6YUOtG6bt26iouLU4sWLWyqyjN69+5d6nVsqlatqjfffFMPPvigl6oCANjNUSM3ixcv1rlz5xQVFaUrr7xSknTgwAF98803+u677zRhwgR16NDB5ioBwNmqVKmiX/3qV+revbv+/Oc/F7qIqgneeOMNLV26VKtWrdK+ffuUnp6urKwshYeHq2XLlurdu7dGjRqlRo0a2V0qAMCLHDVyk5SUpObNmxe6svNXX32lOXPmqE6dOpo9e3aFL8YGAAAAwFyOSgnR0dGFgo106VoUkZGROnHihA4fPmxDZQAAAACczlHhpiR5ozVBQY6aSQcAAADAIXwi3Kxbt06HDx9Wo0aN/O46DQAAAADKxpHDIJ9//rkOHDig7OxsHTp0SAcOHFCdOnU0duxYzrcBAAAAUCRHhputW7dq+/btrt/r1aunMWPGlHnFn/Hjxxe6rUaNGpo6darbagQAAADgLI4MN88884wk6cyZM0pNTdXChQs1ZcoUDR06VHfffbfN1QEAAABwIkctBV2c3NxcTZ48Wfv379df/vIXtWzZslLtHT16VPXq1ZMkZWRkuKNElKJ+/fqS6G9vos+9i/72Lvrb++hz76K/vY8+9668/g4ODnZru44cublcUFCQbrzxRu3bt0+bN2+udLjJn+d8INsZhf72Pvrcu+hv76K/vY8+9y762/voc9/mM2fn16xZU5KUmZlpcyUAAAAAnMhnws3OnTslSQ0bNrS5EgAAAABO5JhpaT/99JNOnDihzp07KzAw0HV7bm6uVq5cqW+++UYhISG68cYbbawSANBt4fwSt68fNMxLlQAAUJBjwk1aWppmzZqlmjVrqnnz5qpZs6aysrKUmpqqEydOKDg4WKNGjXKdfAQA8IzSwktl7k/wAQB4kmPCTbt27XTXXXdp586dSk1NVWZmpoKCghQREaHOnTurX79+ioyMtLtMADBOZcNMZR6LsAMAcCfHhJuIiAgNG8Z/cgDgDd4MNCXJXwdBBwBQWY4JNwAAz3JKoCkOQQcAUFmEGwAwnNNDTVHyaibkAADKg3ADAIbyxVBzOUIOAKA8fOY6NwCAsum2cL4RwSY/054PAMAzGLkBAEOYHgAYxQEAlIaRGwDwcSaO1JTEn54rAKB8GLkBAB/lzwf5bWbPlMQoDgCgIEZuAMAH+XOwyY9+AADkR7gBAB/ib1PQyoL+AADkIdwAgI/gIL549A0AQOKcGwBwPA7cy4bV1AAAjNwAgIMRbMqPPgMA/0W4AQCH4iC94ug7APBPTEsDAIfhwBwAgIph5AYAHIRg4z70JQD4H8INADgEB+PuR58CgH8h3ACAA3AQ7jn0LQD4D8INANiMg2/Po48BwD8QbgDARhx0AwDgPoQbALBBt4XzCTZeRn8DgPkINwDgZRxk24e+BwCzEW4AwIs4uAYAwHMINwDgJQQbZ+B1AABzEW4AwAs4oAYAwPMINwDgYQQb5+E1AQAzEW4AwIM4iAYAwHsINwDgIQQbZ+P1AQDzEG4AwAM4cAYAwPsINwDgZgQb38FrBQBmIdwAgBtxsAwAgH0INwDgJm1mz7S7BFQAgRQAzEG4AQA3INgAAGA/wg0AVBLf/AMA4AyEGwCA3yOgAoAZCDcAUAkcFAMA4ByEGwCoIIINAADOQrgBgAog2AAA4DyEGwAoJ4KNmXhdAcD3EW4AoBw4AAYAwLkINwBQRgQbAACcjXADAGVAsAEAwPkINwBQCoKN/+C1BgDfRrgBgBJwsAsAgO8g3ABAMQg2AAD4FsINABSBYAMAgO8h3ADAZQg2AAD4piC7CwAApyDUAADg2xi5AQARbAAAMAHhBoDfI9gAAGAGwg0Av0awweV4TwCA7yLcAPBbHMQCAGAWwg0Av0SwAQDAPIQbAH6HYAMAgJlYChqA3yDUAABgNkZuAPgFgg0AAOYj3AAwHsEGAAD/QLgBYDSCDQAA/oNzbgAYiVADAID/YeQGgHEINgAA+CfCDQCjEGwAAPBfTEsDYARCDQAAYOQGgM8j2AAAAImRGwA+jFADT1g/aJjdJQAAKoiRGwA+iWADAAAux8gNAJ9CqAEAAMVh5AaAzyDYAACAkjByA8DxCDUAAKAsCDcAHItQAwAAyoNpaQAciWADAADKi5EbAI5CqAEAABVFuAHgCIQaAABQWUxLA2A7gg2cggt4AoBvY+QGgG0INQAAwJ0INwC8jlADAAA8gXADwGsINQAAwJMINwA8jlADAAC8wTHhJjs7W1u3btXmzZu1d+9epaen6+LFi4qMjFTnzp01YMAAVa1a1e4yAZQTwQa+gsUEAMD3OSbcrF+/Xm+99ZYkqUmTJoqJidHZs2f1n//8RwsWLNCGDRs0ZcoU1a5d2+ZKAZQFoQYAAHibY8JNUFCQ+vbtq/79+6tRo0au20+cOKHp06dr//79eu+99zR27FgbqwRQGkINAACwi2PCTY8ePdSjR49Ct9epU0ePPPKIJk+erLi4OOXm5iooyDFlA/gvQg0AALCbT6SEqKgoSdL58+eVlZWlOnXq2FwRgDyEGgAA4BRV7C6gLNLS0iRJgYGBCgsLs7kaAHkINjAFiwkAgBl8YuRm2bJlkqRrr71WwcHBpe4/fvz4QreFhIRo+vTpkqT69eu7prY1aNDAjZWiOPS393myz9vMnun2NgE78dlUOj7HvYv+9j763Ls8dZqJ48NNQkKC1qxZo8DAQA0ZMsTucgC/RqgBAABO5uhwc/DgQb322muyLEv333+/mjZtWqb7zZgxo8TtGRkZql+/viQpPT29smWiDPK+BaG/vcedfc70M5iOz6bS8TnuXfS399Hn3pXX32WZlVUejj3n5vjx43rhhRd05swZDRgwQP369bO7JMAvEWxgOs63AQBzOHLkJjMzU9OmTVNGRoZiY2N1//33210S4HcINQAAwNc4LtycPXtWf/3rX3Xo0CFdf/31GjlypAICAuwuC/AbhBoAAOCrHBVuzp8/r7/97W/au3evYmJiNG7cOFWp4tiZc4BRCDUAAMDXOSY5XLx4Uf/4xz+0Y8cOtW3bVk888YTHlogDUBDBBv6K820AwCyOSQ/Lly9XXFycJKlmzZqaM2dOkfvdf//9qlWrljdLA4xFqAEAACZxTLg5ffq06+e8kFOUe++9l3ADVBKhBgAAmMgx4Wbw4MEaPHiw3WUAxiPYAAAAUzkm3ADwLEINUBDn2wCAeRyzoAAAzyHYAAAAf8DIDWCwNrNn2l0CAACA1zByAxiKYAMUjylpAGAmRm4AwzAFDQAA+CtGbgCDEGwAAIA/Y+QGMAChBgAAgJEbwOcRbIDy4XwbADAX4QbwYQQbAACA/2FaGuCDCDUAAACFMXID+BiCDVBxTEkDALMRbgAfQrABAAAoHtPSAB9AqAEAACgdIzeAwxFsAPdgShoAmI9wAzgYwQYAAKDsCDeAQxFsAAAAyodwAzgQwQYAAKD8WFAAcBBCDeAZnG8DAP6BkRvAIQg2AAAAlUO4ARyAYAMAAFB5hBvAZgQbwLOYkgYA/oNwA9iIYAMAAOA+hBvAJgQbAAAA9yLcADYg2AAAALgf4QbwMoIN4D2cbwMA/oVwA3gRwQYAAMBzCDeAlxBsAAAAPItwA3gBwQYAAMDzCDeAhxFsAHtwvg0A+B/CDeBBBBsAAADvIdwAHkKwAQAA8C7CDeABBBsAAADvI9wAbkawAezH+TYA4J8IN4AbEWwAAADsQ7gBAAAAYATCDeAmjNoAAADYi3ADuAHBBgAAwH6EG6CSCDaAs7CYAAD4L8INUAkEGwAAAOcg3AAAAAAwAuEGqCBGbQAAAJyFcANUAMEGAADAeQg3QDkRbADnYjEBAPBvhBsAAAAARiDcAOXAqA0AAIBzEW6AMiLYAAAAOBvhBgAAAIARCDdAGTBqAzgfiwkAAAg3QCkINgAAAL6BcAMAAADACIQboASM2gAAAPgOwg1QDIINAACAbyHcAAAAADAC4QYoAqM2AAAAvodwAwAAAMAIhBvgMozaAL6Ha9wAACTCDVAAwQYAAMB3EW4AAAAAGIFwA/wXozYAAAC+jXADAAAAwAiEG0CM2gAAAJiAcAMAAADACIQb+D1GbQAAAMxAuAEAAABgBMIN/BqjNoDv4wKeAIA8hBsAAAAARiDcwG8xagMAAGAWwg0AAAAAIxBuAAAAABiBcAO/xJQ0AAAA8xBuAAAAABiBcAO/w6gNAACAmQg3AAAAAIxAuIFfYdQGAADAXIQbAAAAAEYIsruA/Pbt26dt27Zpz5492r17t06cOKHg4GB9+OGHdpcGAAAAwOEcFW4WLlyoH374we4yYCimpAEAAJjNUeGmdevWatq0qVq0aKEWLVro0UcftbskAAAAAD7CUeFm4MCBdpcAAAAAwEexoAD8AlPSAAAAzEe4AQD4rPWDhtldAgDAQQg3AAAAAIzgqHNu3GX8+PGFbgsJCdH06dMlSfXr11dQ0KWn3qBBA6/W5q/s7O82s2d6/TEBeAef4d7D/5veRX97H33uXXn97W6M3AAAAAAwgpEjNzNmzChxe0ZGhurXry9JSk9P90ZJfi/vWxD6G4A78ZniPXyOexf97X30uXfl9XdwcLBb22XkBkZjlTQAAAD/QbgBAAAAYATCDQAAAAAjEG5gLKakAQAA+BdHLSiQkJCgRYsWFbgtNzdXkyZNcv1+zz336LrrrvN2aQAAAAAczlHhJjMzU7t37y5wm2VZBW7LzMz0dlkAAAAAfICjwk1sbKxiY2PtLgMAAACAD+KcGxiJ820AAAD8D+EGAAAAgBEINwAAAACMQLgBAAAAYATCDYzD+TYAAAD+iXADAAAAwAiEGwAAAABGINwAAAAAMALhBkbhfBsAAAD/RbgBAAAAYATCDQAAAAAjEG4AAAAAGIFwAwAAAMAIhBsYg8UEAAAA/BvhBgAAAIARCDcAAAAAjEC4AQAAAGAEwg0AAAAAIxBuAAAAABiBcAMjsFIaAAAACDcAAAAAjEC4AQAAAGAEwg0AAAAAIxBuAAAAABiBcAMAAADACIQb+DxWSgMAAIBEuAEAAABgCMINAAAAACMQbgAAAAAYgXADAAAAwAiEGwAAAABGINwAAAAAMALhBgAAAIARCDcAAAAAjEC4gU/jAp4AAADIQ7gBAAAAYATCDQAAAAAjEG4AAAAAGIFwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEG/isbgvn210CAAAAHIRwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEGwAAAABGINwAAAAAMALhBgAAAIARCDcAAAAAjEC4AQAAAGAEwg0AAAAAIxBuAAAAABiBcAMAAADACIQbAAAAAEYg3AAAAAAwAuEGAAAAgBEINwAAAACMQLgBAAAAYATCDQAAAAAjEG4AAAAAGIFwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEYLsLuByOTk5WrJkiTZs2KCMjAyFhYUpJiZGQ4YMUb169ewuDwAAAIBDOWrkJicnR1OnTtXChQt17tw5dezYUfXq1dPatWv15JNP6ujRo3aXCAAAAMChHDVys3jxYu3atUutW7fW5MmTVbVqVUnS0qVL9f7772v27Nl6/vnnba4STtBt4Xy7SwAAAIDDOGbkJjc3V8uXL5ckPfLII65gI0kDBgxQVFSUfvrpJ+3bt8+uEgEAAAA4mGPCTVJSks6cOaOGDRuqWbNmhbZ37txZkvTDDz94uzQAAAAAPsAx4SYlJUWSigw2ktS8efMC+wEAAABAfo455yYjI0OSil0RrW7dugX2K8n48eML3RYSEqLp06dLkurXr6+goEtPvUGDBhWqF+VDfwPwBD5TvIfPce+iv72PPveuvP52N8eM3Jw7d06SFBoaWuT2vHNw8vYDAAAAgPwcM3JjWValtuc3Y8aMErdnZGSofv36kqT09PQyt4uKy/sWhP4G4E58pngPn+PeRX97H33uXXn9HRwc7NZ2HTNyU61aNUlSdnZ2kdvzbs+/ihoAAAAA5HFMuMkbSTl+/HiR23/++ecC+wEAAABAfo4JN1FRUZKk/fv3F7k97/o2efsBAAAAQH6OCTfR0dGqXr260tLSigw4mzZtkiRdd9113i4NDrR+0DC7SwAAAIDDOCbcBAUF6dZbb5UkzZ07t8CqaEuXLlVKSoqio6PVsmVLu0oEAAAA4GCOWS1Nku6++25t375du3bt0tixYxUdHa2MjAzt3r1bNWvW1KhRo+wuEQAAAIBDOWbkRrp0oc3nnntO99xzj0JCQhQfH69jx46pR48eevHFFxUZGWl3iQAAAAAcylEjN9KlgDNkyBANGTLE7lIAAAAA+BBHjdwAAAAAQEURbgAAAAAYgXADAAAAwAiEGwAAAABGINwAAAAAMALhBgAAAIARCDcAAAAAjEC4AQAAAGAEwg0AAAAAIxBuAAAAABiBcAMAAADACIQbAAAAAEYg3AAAAAAwAuEGAAAAgBEINwAAAACMQLgBAAAAYATCDQAAAAAjEG4AAAAAGIFwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEG/is9YOG2V0CAAAAHIRwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEGwAAAABGINwAAAAAMALhBj5t/aBhdpcAAAAAhyDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEGwAAAABGINzA53EhTwAAAEiEGwAAAACGINwAAAAAMALhBgAAAIARCDcAAAAAjEC4AQAAAGAEwg0AAAAAIxBuYASWgwYAAADhBgAAAIARCDcAAAAAjEC4AQAAAGAEwg0AAAAAIxBuYAwWFQAAAPBvhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuIFRWFQAAADAfxFuAAAAABiBcAMAAADACIQbAAAAAEYg3AAAAAAwAuEGxmFRAQAAAP9EuAEAAABgBMINAAAAACMQbgAAAAAYgXADI3HeDQAAgP8JsrsASTp37pzi4uK0Z88e7d69WykpKcrNzdXw4cM1cOBAu8sDAAAA4AMcEW6OHj2q119/3e4yAAAAAPgwR4SbqlWrqmfPnmrZsqVatGihTZs26dNPP7W7LPi49YOGqdvC+XaXAQAAAC9xRLiJjIzUyJEjXb/Hx8fbWA0AAAAAX8SCAgAAAACMQLgBAAAAYATCDYzGktAAAAD+wxHn3Ljb+PHjC90WEhKi6dOnS5Lq16+voKBLT71BgwZerc1f0d8APIHPFO/hc9y76G/vo8+9K6+/3d6uOxp55ZVXdODAgXLdZ/To0WrZsqU7Hh4AAAAA3BNu0tPTdfjw4XLdJzs72x0PXaQZM2aUuD0jI0P169eXdKl2eF7etyB29DdLQgPm4jPce+z8HPdH9Lf30efeldffwcHBbm3XLeEmb7oXAADe1G3hfM6tAwC4sKAAAAAAACMQbuAX+GYXAADAfIQbAAAAAEZwzFLQL730kk6ePClJOn78uCRpxYoVio+PlySFh4drwoQJdpUHA7CwAAAAgNkcE26Sk5MLrU5x/PhxV9BhzXEAAAAAJXFMuHnjjTfsLgEAAACAD+OcG/gVFhYAAAAwF+EGAAAAgBEIN/A7jN4AAACYiXADAAAAwAiEGwAAAABGINzALzE1DQAAwDyEGwAAAABGINzAbzF6AwAAYBbCDQAAAAAjEG7g1xi9AXxft4Xz7S4BAOAQhBsAAAAARiDcwO8xegMAAGAGwg0AAAAAIxBuADF6AwAAYALCDQAAAAAjEG6A/2L0BgAAwLcRbgAAAAAYgXAD5MPoDQAAgO8i3AAAAAAwAuEGuAyjN4Dv6bZwvt0lAAAcgHADFIGAAwAA4HsINwAAAACMQLgBisHoDQAAgG8h3AAAAAAwAuEGKAGjNwAAAL6DcAOUgoADAADgGwg3AAAjsBw0AIBwA5QBozcAAADOR7gByoiAAwAA4GyEGwAAAABGINwA5cDoDQAAgHMRboByIuAAAAA4E+EGqAACDuBMrJgGAP6NcAMAAADACIQboIIYvQEAAHAWwg1QCQQcAAAA5yDcAJVEwAEAAHAGwg0AwCgsKgAA/otwA7gBozcAAAD2I9wAbkLAAQAAsBfhBnAjAg4AAIB9CDeAmxFwAAAA7EG4ATyAgAPYi0UFAMA/EW4ADyHgAAAAeBfhBvAgAg4AAID3EG4AAAAAGIFwA3gYozeAPTjvBgD8D+EG8AICDgAAgOcRbgAvIeAAAAB4FuEG8CICDuBdTE0DAP9CuAG8jIADAADgGYQbwAYEHAAAAPcj3AA2WT9oGCEHAADAjQg3gM0IOIBncd4NAPgPwg3gAAQcAACAyiPcAA5BwAEAAKgcwg3gIAQcAACAiiPcAA5DwAHcj/NuAMA/EG4AByLgAAAAlB/hBnAolooGAAAoH8IN4HAEHMA9mJoGAOYj3AA+gIADAABQOsIN4CMIOAAAACUj3AA+hPNwAAAAike4AXwQAQeoGM67AQCzEW4AH0XAAQAAKIhwA/gwpqkBAAD8D+EGMAABByg7pqYBgLkIN4AhGMUBAAD+jnADGIaAAwAA/BXhBjDQ+kHDtOsPj9ldBuBYTE0DADMRbgCDEXAAAIA/IdwAhuNcHAAA4C8IN4CfIOAAAADTEW4AP8IoDvA/nHcDAOYJsrsASTp06JDi4+O1bds2HTlyRKdOnVKNGjXUpk0b9e/fX23btrW7RMAo6wcN48AOAAAYxxHhZurUqfr5559VrVo1tWrVSq1bt9bBgwcVFxen+Ph4PfDAA+rfv7/dZQJGyRvBIeQAAABTOCLcNG7cWPfdd59uuOEGBQX9r6SVK1fqnXfe0QcffKCYmBg1btzYxioBMxFy4M+6LZzPVE0AMIgjzrmZPHmyunXrViDYSFKfPn0UExOjixcvauPGjTZVB/gHDvAAAICvc0S4KUlUVJQk6cSJEzZXApiPBQcAAIAvc3y4SUtLkySFh4fbWwjgRwg5AADAFzk63Bw9elQJCQmSpI4dO9pcDeB/CDnwB5xvBgDmCLAsy7K7iKJcuHBBzz//vJKSknTjjTdq3LhxZb7v+PHjC90WEhKi6dOnS5LOnz/vOr8nNzfXLfWiZPS397m7z9vMnumWdgAn2vWHx+wuwfH4HPcu+tv76HPvyuvvgIAA97brjkZeeeUVHThwoFz3GT16tFq2bFns9rlz5yopKUkNGzbUiBEjKlsigErKO/gj5AAAAKdyS7hJT0/X4cOHy3Wf7OzsYrctXLhQK1euVO3atTVp0iSFhYWVq+0ZM2aUuD0jI0P169eXdKl2eF6DBg0k0d/e5Kk+Z+lomIjPptLxOe5d9Lf30efeldffwcHBbm3XLeEmb7qXOyxfvlwLFixQ9erVNWnSJEVGRrqtbQDuQ8iBSbjeDQCYwVELCnz77beaN2+eQkNDNXHiRDVt2tTukgCUggNCAADgFI4JNwkJCZo1a5YCAwP1xBNPKDo62u6SAJQRq6oBAAAncES4SUpKcp0nM27cOMXExNhcEYCKIOQAAAA7ueWcm8p68cUXlZOTo4iICMXHxys+Pr7QPtHR0erVq5cN1QEoL87HAQAAdnBEuDlz5owk6dixYzp27Fix+xFuAN9CyIEvYVEBAPB9jgg3CxYssLsEAB60ftAwAg4AAPA4R5xzA8B8nI8DAAA8jXADwKsIOQAAwFMINwBsQcgBAADuRrgBYCtCDgAAcBfCDQBHIODACVj4AgB8G+EGgGMwigMAACqDcAPAcQg5AACgIgg3AByLgAMAAMqDcAPA0RjFAQAAZUW4AeATCDgAAKA0hBsAPoNRHAAAUBLCDQCfQ8gBAABFIdwA8FkEHHgC17oBAN9FuAHg0wg4AAAgD+EGgM9jmhoAAJAINwAMQsABAMC/EW4AGIVRHAAA/BfhBoCRCDgAAPgfwg0AYxFwAADwL4QbAEZjmhoAAP6DcAPALxBwAAAwH+EGgN8g4AAAYDbCDQC/QsABAMBchBsAfoeAAwCAmQg3APwSCw0AAGAewg0Av0bAAQDAHIQbAH6PgIP8eD8AgO8i3ACAOKAFAMAEhBsA+C8CDgAAvo1wAwD5EHAAAPBdhBsAuAwBBwAA30S4AYAiEHAAAPA9hBsAKAYBBwAA30K4AYASEHD8C683APg2wg0AlIIDXgAAfAPhBgDKgIADAIDzEW4AoIwIOAAAOBvhBgDKgYADAIBzEW4AoJwIOGbidQUA30e4AYAK4EAYAADnIdwAQAURcAAAcBbCDQBUAgHHDLyOAGAGwg0AVBIHxgAAOAPhBgDcYNcfHrO7BAAA/B7hBgDchIDjmxh5AwBzEG4AwI04UAYAwD6EGwBwMwKO7+C1AgCzEG4AwAM4aAYAwPsINwDgIQQcAAC8i3ADAB5EwHEuXhsAMA/hBgA8jINoAAC8g3ADAF5AwHEWXg8AMBPhBgC8hANqZ+B1AABzEW4AwIs4sAYAwHMINwDgZQQc+9D3AGA2wg0A2ICDbAAA3I9wAwA2IeB4F/0NAOYj3ACAjTjg9g76GQD8A+EGAGzGgbdn0b8A4D8INwDgAOsHDeMgHACASiLcAICDEHDci/4EAP9CuAEAh+GA3D3oRwDwP4QbAHAgDswrh/4DAP9EuAEAh+I8nIqhzwDAfxFuAMDhOFgvO/oKAPwb4QYAfACjOKWjfwAAhBsA8CEcwBeNfgEASFKQ3QUAAMon70C+28L5NldiP0INACA/Rm4AwEf5+4H9rj88ZncJAACHIdwAgA/z13NxCDYAgKIwLQ0ADOAvU9X8McgBAMqOcAMABjE15BBqAABlwbQ0ADCQSdPVTHkeAADPY+QGAAzmyyM5hBoAQHkRbgDAD+QPCk4OOgQaAEBlEG4AwM84LegQaAAA7uKIcJOSkqJVq1Zp3759ysjIUFZWloKDg9W4cWPddNNN6tOnjwIDA+0uEwCMY1fQIdAAADzBEeHmp59+0ooVK9SgQQM1btxYNWvWVFZWlpKSkrR7927FxcXp6aefVlCQI8oFACMVFzgqG3oIMgAAb3FEWujQoYNee+01NWzYsMDtJ0+e1NSpU/Xjjz9q9erV6tu3r00VAoD/IpwAAHyFI5aCbtiwYaFgI0nh4eEaOHCgJOnHH3/0clUAAAAAfIkjwk1JqlS5VCJT0gAAAACUxNGJ4fTp01q6dKmkS1PX3CUgIKDIn+F59Lf30efeRX97F/3tffS5d9Hf3kef+7YAy7Isu4vIc+TIEX366aeyLEunTp3Srl27dO7cOfXu3Vu/+93vyvxmGz9+fKHbatSooalTp7q7ZAAAAAAO4aiRm1OnTmndunUFbrv11ls1dOhQUjQAAACAErll5OaVV17RgQMHynWf0aNHq2XLlkVuu3jxojIyMhQXF6dPPvlE4eHhmjRpkiIiIipbqsvEiRMlSdOnT3dbmyge/e199Ll30d/eRX97H33uXfS399Hn3uWp/nbLyE16eroOHz5crvtkZ2cXu61KlSqKiIjQgAEDFBERoZdffllz5851dYI75OTkuK0tlI7+9j763Lvob++iv72PPvcu+tv76HPv8lR/uyXceDLhdurUSVWrVlViYqJyc3NZNQ0AAABAkRy/FHRAQIDCwsJ08eJFnT592u5yAAAAADiU48NNWlqajh8/rmrVqqlWrVp2lwMAAADAoRwRbj777DOlpaUVuv3w4cOaOXOmLMtSjx49XBf0BAAAAIDLOeI6N3/84x+VkZGhpk2bqmHDhpIuLVKwb98+WZaltm3b6qmnnlLVqlVtrhQAAACAUzki3Hz77bfasmWL9u7dq5MnTyonJ0dhYWFq2rSpunbtqu7duzNqAwAAAKBEjgg3AAAAAFBZDIcAAAAAMALhBgAAAIARCDcAAAAAjEC4AQAAAGAEwg0AAAAAIwTZXYCTLVy4UAsWLJAkjR07Vl27drW5InOkpKRo1apV2rdvnzIyMpSVlaXg4GA1btxYN910k/r06aPAwEC7yzTGoUOHFB8fr23btunIkSM6deqUatSooTZt2qh///5q27at3SUa59y5c4qLi9OePXu0e/dupaSkKDc3V8OHD9fAgQPtLs9n5eTkaMmSJdqwYYMyMjIUFhammJgYDRkyRPXq1bO7PKPs27dP27Ztc72HT5w4oeDgYH344Yd2l2ak7Oxsbd26VZs3b9bevXuVnp6uixcvKjIyUp07d9aAAQO43p8HLF26VElJSUpNTdWpU6d0/vx5hYeHq127drrzzjvVpEkTu0s01unTpzVu3DhlZmbqiiuu0KuvvuqWdgk3xTh8+LAWL16sgIAAsVq2+/30009asWKFGjRooMaNG6tmzZrKyspSUlKSdu/erbi4OD399NMKCuIt6g5Tp07Vzz//rGrVqqlVq1Zq3bq1Dh48qLi4OMXHx+uBBx5Q//797S7TKEePHtXrr79udxlGycnJ0dSpU7Vr1y7VqVNHHTt2VHp6utauXauEhARNmzZNkZGRdpdpjIULF+qHH36wuwy/sX79er311luSpCZNmigmJkZnz57Vf/7zHy1YsEAbNmzQlClTVLt2bZsrNcvixYt17tw5RUVF6corr5QkHThwQN98842+++47TZgwQR06dLC5SjP985//VFZWltvb5cixCJZl6a233lL16tXVqlUrPtw9oEOHDnrttdfUsGHDArefPHlSU6dO1Y8//qjVq1erb9++NlVolsaNG+u+++7TDTfcUCAwrly5Uu+8844++OADxcTEqHHjxjZWaZaqVauqZ8+eatmypVq0aKFNmzbp008/tbssn7Z48WLt2rVLrVu31uTJk13fYi9dulTvv/++Zs+ereeff97mKs3RunVrNW3aVC1atFCLFi306KOP2l2S0YKCgtS3b1/1799fjRo1ct1+4sQJTZ8+Xfv379d7772nsWPH2lileSZMmKDmzZsrJCSkwO1fffWV5syZozfffFOzZ8/mYvJutn37dq1bt069e/fWqlWr3No2r1QRvv76a/3000964IEHVKNGDbvLMVLDhg0LBRtJCg8Pd03Z+fHHH71clbkmT56sbt26FRoJ69Onj2JiYnTx4kVt3LjRpurMFBkZqZEjR6p3795q1qwZ/zFWUm5urpYvXy5JeuSRRwpMzxkwYICioqL0008/ad++fXaVaJyBAwdq8ODB+vWvf63w8HC7yzFejx49NGLEiALBRpLq1KmjRx55RJIUFxen3NxcO8ozVnR0dKFgI0l9+/ZVZGSkTpw4ocOHD9tQmblycnL0zjvvqHHjxrr99tvd3j7/217m5MmT+vDDD3X11Vfrpptusrscv5R3EMiUNO+IioqSdOnbQcCpkpKSdObMGTVs2FDNmjUrtL1z586SxEg7jJT3OX3+/HmPTONB0Tge8YxPPvlEaWlpGjFihEfOrybcXGbu3LnKycnRiBEj7C7FL50+fVpLly6VJOa4eklaWpok8c0sHC0lJUWSigw2ktS8efMC+wEmyfucDgwMVFhYmM3V+Id169bp8OHDatSokSIiIuwuxxgpKSlaunSpYmNj1a5dO488BlE0n82bN+v777/X4MGDCw0LwzOOHDmiTz/9VJZl6dSpU9q1a5fOnTun3r17q1u3bnaXZ7yjR48qISFBktSxY0ebqwGKl5GRIUnFrohWt27dAvsBJlm2bJkk6dprr1VwcLDN1Zjp888/14EDB5Sdna1Dhw7pwIEDqlOnjsaOHcu0Yje5ePGi65z2++67z2OPQ7j5r3PnzmnOnDlq1KiR7rzzTrvL8RunTp3SunXrCtx26623aujQoQoICLCpKv9w4cIFzZo1S+fPn9eNN97o+uYbcKJz585JkkJDQ4vcnncOTt5+gCkSEhK0Zs0aBQYGasiQIXaXY6ytW7dq+/btrt/r1aunMWPG8H+jGy1fvlx79uzRqFGjVLNmTY89jjHh5pVXXtGBAwfKdZ/Ro0erZcuWkqSPPvpIx48f17PPPsu3ImVQ2f7OEx0drQULFujixYvKyMhQXFycPvnkE23dulWTJk1iKPi/3NXf+c2dO1dJSUlq2LAh0zCL4Ik+R8WVtiQ/S/bDRAcPHtRrr70my7J0//33q2nTpnaXZKxnnnlGknTmzBmlpqZq4cKFmjJlioYOHaq7777b5up8X0ZGhj7++GO1a9dOsbGxHn0sY8JNenp6uVezyM7OliTt2bNHK1asUPfu3XXVVVd5ojzjVKa/i1KlShVFRERowIABioiI0Msvv6y5c+dq4sSJlS3VCO7u74ULF2rlypWqXbu2Jk2axBzuIri7z1E51apVk1R8H+fdzkUOYYrjx4/rhRde0JkzZzRgwAD169fP7pL8Qo0aNdS2bVs99dRTmjx5sv71r3/pmmuu4YurSpozZ45yc3O98mWqMeFm+vTpFb5vQkKCLMtSamqqpkyZUmDboUOHJP3vYPCGG27QrbfeWplSjVCZ/i5Np06dVLVqVSUmJio3N5dVSuTe/l6+fLkWLFig6tWra9KkSVz0sBiefI+j/OrXry/p0gFfUX7++ecC+wG+LDMzU9OmTVNGRoZiY2N1//33212S3wkKCtKNN96offv2afPmzYSbSkpISFCNGjU0Z86cArefP39e0qWRnbxj8IkTJ1bqiyqOGvNJTk4udtuhQ4d06NAhhoS9ICAgQGFhYcrIyNDp06dZxcuNvv32W82bN0+hoaGaOHEi72f4jLylcPfv31/k9rzr2+TtB/iqs2fP6q9//asOHTqk66+/XiNHjuQcVJvknReSmZlpcyVmOHPmjHbu3FnktpycHNe2CxcuVOpxCDeSBg8erMGDBxe57Y033tC6des0duxYde3a1cuV+ae0tDQdP35c1apVU61atewuxxgJCQmaNWuWAgMD9cQTTyg6OtrukoAyi46OVvXq1ZWWlqb9+/cXWhJ606ZNkqTrrrvOjvIAtzh//rz+9re/ae/evYqJidG4ceNYqctGeQfbRV10HOWzYMGCIm8/duyYRo8erSuuuEKvvvqqWx6LvxjY4rPPPnOt25/f4cOHNXPmTFmWpR49evCh7iZJSUmaMWOGJGncuHGKiYmxuSKgfIKCglxTgufOnVtgVbSlS5cqJSVF0dHRTB2Bz7p48aL+8Y9/aMeOHWrbtq2eeOIJpmV72E8//aTvvvuu0EhBbm6uvvzyS33zzTcKCQnRjTfeaFOFqAj+amCLr776Sh999JGaNm3q+kYkPT1d+/btk2VZatu2rYYPH25zleZ48cUXlZOTo4iICMXHxys+Pr7QPtHR0erVq5cN1ZnrpZde0smTJyX971yRFStWuPo/PDxcEyZMsKs8n3P33Xdr+/bt2rVrl8aOHavo6GhlZGRo9+7dqlmzpkaNGmV3iUZJSEjQokWLCtyWm5urSZMmuX6/5557GC1zk+XLlysuLk7SpelQl5+bkOf+++9nVoObpKWladasWapZs6aaN2+umjVrKisrS6mpqTpx4oSCg4M1atQozuXzMYQb2GLo0KHasmWL9u7dq61btyonJ0dhYWG65ppr1LVrV3Xv3p1RGzc6c+aMpEvDv8eOHSt2P8KNeyUnJys9Pb3AbcePH3cFnQYNGthRls8KCQnRc889p8WLF2v9+vWKj49XjRo11KNHDw0ZMoQDEDfLzMzU7t27C9xmWVaB2zgXwX1Onz7t+jkv5BTl3nvvJdy4Sbt27XTXXXdp586dSk1NVWZmpoKCghQREaHOnTurX79+LLrjgwIsLg4AAAAAwAB8NQ4AAADACIQbAAAAAEYg3AAAAAAwAuEGAAAAgBEINwAAAACMQLgBAAAAYATCDQAAAAAjEG4AAAAAGIFwAwAAAMAIhBsAAAAARiDcAAAAADAC4QYAAACAEQg3AAAAAIxAuAEAAABgBMINAAAAACMQbgAAAAAYgXADAAAAwAiEGwAAAABG+P/bGnMQ+4/2WQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Coefficient for classical RK4\n",
    "c2, c3, c4 = 0.5, 0.5, 1.0\n",
    "a21 = 0.5\n",
    "a32 = 0.5\n",
    "a43 = 1.0\n",
    "b1, b2, b3, b4 = 1/6, 1/3, 1/3, 1/6\n",
    "\n",
    "k1 = ff(t, y)\n",
    "k2 = ff(t + c2*h, y+ h*a21*k1)\n",
    "k3 = ff(t + c3*h, y+ h*a32*k2)\n",
    "k4 = ff(t + c4*h, y+ h*a43*k3)\n",
    "\n",
    "# Factor expression\n",
    "expr = sy.simplify((y+ h*(b1*k1 + b2*k2 + b3*k3 + b4*k4)) / y)\n",
    "print(expr)\n",
    "\n",
    "# Sig\n",
    "z = sy.Symbol('z')\n",
    "sigf = sy.lambdify(z, expr.subs(h*lam, z), 'numpy')\n",
    "\n",
    "# Make grid\n",
    "xx = np.linspace(-3, 3, 201)\n",
    "yy = np.linspace(-3, 3, 201)\n",
    "X, Y = np.meshgrid(xx, yy)\n",
    "z = X + Y*1j\n",
    "\n",
    "# Amplication factor of Explicit Euler (z=lambda h)\n",
    "sig = sigf(z)\n",
    "\n",
    "# Same scale of x and y\n",
    "plt.axis('equal')\n",
    "\n",
    "# Stability region\n",
    "plt.contourf(X,Y,abs(sig), levels=[0, 1])\n",
    "\n",
    "# Title\n",
    "plt.title('Stability region of RK4')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}